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ABSTRACT 



In this thesis we deal with the generation of bivariate 
random variables with identical marginal distributions and 
specified correlation coefficients. Many schemes have been 
put forward for different cases; they almost always involve 
computing the inverse probability distribution of the 
marginal distributions or exploiting special properties of 
the random variables in question. 

A simple mixture-truncation method is put forward here 
to generate bivariate distributions with any marginal 
distribution from univariate generators. 
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I . INTRODUCTION 



In many simulation applications, it is required to 
generate bivariate random variables which have identical 
marginal distribution. 

For example, in reliability problems, the assumption that 
two components have independent exponential failure times is 
very often unrealistic, and much effort has gone into deriving 
bivariate exponential random variables to handle these situa- 
tions [Gaver (1972) , Olkin & Marshall (1967) , Downton (1970) ] . 
Now except in very specific physical situations it may be 
difficult to specify the complete bivariate distribution of 
life time of each component. However, it may be realistic 
to specify the marginal distributions and some measure of 
dependence (usually the correlation coefficient) between the 
life time of each component. In this kind of situation, we 
can use bivariate random vectors having given marginal dis- 
tribution and dependence to solve the problem in simulation. 

To generate these vectors, there exists some previous works 
and existing methods, discussed in Section II. As seen in 
Section II, most previous work is specific to specified 
marginal distributions and uses inverse transformation methods 
as a basic concept. 

An example is the recent work by Johsnon and Tenenbein 
(1979), to generate a bivariate random vector (X,Y) which has 
marginal distribution F^^(x), F 2 (y) and correlation p , by a 
weighted linear combination method. 
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Define 



X = F^^(H^(U)) 

Y = F2^(H2(V)) 



or 



Y = F^^d - (V) ) 

where and H 2 are the cumulative distribution functions 
(c.d.f.) of U and V respectively and 



U = U' 

V = cU' + (1 - c) V 

where U', V are i.i.d. random variables with probability 
density functino g( ). In this procedure, ^ 

are specific to the marginal distribution and correlation 
desired. The functions F^^ and F 2 ^ are difficult to compute 
in most cases and the weighting factor c is also difficult to 
calculate. Moreover most of the work in univariate random 
number generation has been aimed at avoiding having to calcu- 
late inverse cumulative distribution functions such as 
F^^(.) and F 2 ^(*). These are the reasons why many proposed 
methods are specific to a specified marginal distribution. 
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Again special properties of certain random variables such 
as infinite divisibility have been exploited to give easily 
generated bivariate random variables, often though with 
limited ranges of dependency. One very clever scheme by 
Gaver (1972) to generate bivariate exponential random varia- 
bles uses the fact that the sum of a geometrically distributed 
number of exponential random variables (Y) is exponentially 
distributed and that the minimum of this geometrically dis- 
tributed number of independent logistic random variables (Z) 
is exponential. Clearly when Y is large, Z is small. This 
scheme is of course very specific to exponential marginal 
distributions and, via an exponential transformation, to 
uniform random variables. To avoid these kinds of limitations 
and to make the generation of bivariate random variables 
simpler and more automatic in simulations we develop here a 
scheme presented in Jacobs and Lewis (1977) . 

This scheme, the mixture-truncation method, is a very 
general tool which requires only that a method be available 
for generating random variables with the desired marginal 
distribution. The mixture- truncation method scheme for 
generating bivariate random variables is as follows. 

Let F(x) be the common marginal distribution, of the 
bivariate random variable (Y,Z). Let p be the desired 
correlation between Y and Z (which may or may not be attain- 
able generally or in particular with the mixture-truncation 
scheme) . 
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Define the transition matrix P as 



P = 



1 - a. 






with stationary vector 



1 - a. 



1 - a. 



ir P = ir = ( 



l-aj^ + l~a2^ l-aj^ + l~a2‘ 



and let the range of the random variable X with distribution 
function F(x) be x* Then generate (Y,Z) as follows. 

1 . ( Initiallization) 

i) Choose an "allowable" x^ from (x,^,x^) (the 

allowable range of x^ will usually be smaller 
than X and depends on p and F(x)). 

ii) Set = F(x^), tt 2 = ^ ” ''^1 compute and 

a. 2 • 

iii) Denote by the random variable X truncated to the 

i 

left of x^ and by X^ truncated to the right of 

X . 

o 

2 . (Generate Y) . Choose Y from X^^ with probability 
or choose Y from X2 with probability t ^ 2 ' 

3 . (Generate Z) 

i) If Y is chosen from X^^, choose Z from X^^ with 

probability or choose Z from X2 with probability 
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ii) If Y is chosen from X^, choose Z from with 
probability 1 ~ a2 or choose Z from X2 with 
probability a2* 

In Section III, we will show that the correlation p between 
Y and Z, if is properly chosen, is 



p 


jZ-E(Z), 

°Y °Z 




e M, 



where 



6 = 


- (1 - 02 ) , 


M = 


( ^1 - ^ 2 ^ ^ ^1 ^2 
2 

'^X 


'^1 


= E{X^}, M2 = ' 


2 

^X 


= VAR [X ] , 



Moreover, the generated bivariate random vector (Y, Z) will 
have marginal distributions F (x) and correlation p and the 
joint distribution function of (Y,Z) will be 
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if 



F(z) F(y) 



z < X , Y < X 
— o — o 



(1-ct,) 

+ — [F(z) -F (x^) ] }F (y) if z > Xq/ yi^o 

F(y,z) = ,1. ) 

{a, + [F(y)-F(x^) ] }F(z) if z<x , y>x 

1 Tl 2 ^ — o o 



afTTi + (l-Oi) [F (z) -F(X q) ] — if 2 < x^, y > 



+ {( 1 -a.) + — [F(z)-F(x^) ] } [F(y)-F(x ) ] 
Z TT 2 O o 



These relationships will be developed in Section III. Note 
that (Y,Z) may be continuous or discrete random variables 
or mixtures of both, though in this thesis we concentrate on 
continuous cases. 

The key problem in this very simple algorithm comes at 
the initialization steps i) and ii) . There are two degrees 
of freedom in the selection matrix P but setting = F(x^) 
constraints reduce to one degree of freedom. Specifying 
a desired correlation further constrains the degree of free- 
dom, though not completely. Subject to the constraint that 



and 


a 2 probabilities there may 


be 


a. 


no values x^ which will give 


(y,Z) with correlation 


b. 


P 

one value x 






o 




c . 


a range Xq- 





12 



The main numerical problem of initialization step (i) then 
is to compute the range of allowable for a given p and 
F(x) . 

The main statistical problem is then to choose which of 

the bivariate random vectors (Y,Z) indexed by x^ e Xq^o use. 

An alternate solution to the statistical problem, giving 

another algorithm is then to let x^ have some distribution in 

the range Xq' possibly uniform or triangular, this not only 

alleviates the problem of picking a particular x^ e 

it also smoothes out the distribution of (Y,Z) and possibly 

makes it continuous. The computation of x given p is 

o 

illustrated in subsequent sections for uniform, exponential 
and gamma marginals for (Y,Z). 

Before doing this and developing, in Section III, the 
results already given here, we review a few existing 
methods for generating bivariate random variables in Section 
II. These are methods which seem fairly tractable, for use 
on a computer. 
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II. REVIEW OF EXISTING METHODS 



There are several existing methods for generating 
bivariate random vectors, but most of them are specific to 
particular marginal distribution and use inverse transfor- 
mation method. The inverse transformation method is a very 
useful univariate procedure which, unfortunately, is not 
possible to use with many distributions because it is diffi- 
cult and/or uneconomic to compute the inverse functions. We 
survey here some of the methods which are germane to this 
thesis, in particular concentrating on bivariate random 
variables with uniform, exponential and gamma marginals. First 
we review the problem of determining the range of correlation 
coefficient p which can be obtained for bivariate distribution 
with given marginal distributions, again considering only 
the continuous case. 

A. RANGE OF CORRELATION COEFFICIENT p 

Suppose that Y, Z are random variables with an arbitrary 
joint distribution F(y,z) with finite second moments. Then 
in general the correlation coefficient p can take all values 
in the closed interval [-1,1]. But with specific marginal 
distribution Fj^(y) and F 2 (z), the class of all F(y,z) need 
.not attain the values of -1, 1, of p. The necessary and 
sufficient conditions that there exist determination of F(y,z) 
with p equal to 1 and -1 are given by Moran (1967) as follows. 
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i) there exist constants a and 3 such that aY + 3 has 
the same distribution as Z, 
ii) the distribution of Z is symmetrical about its 
mean. 

To see this, rescale Y to have the same mean and variance 

as Z. If there exists an F(y,z) such that p = 1 we have 

2 . . 

E[Z-Y] = 0, SO that Y = Z with probability one. Then Y 

must have the same distribution as Z. On the other hand 

if there exists an F(y,z) such that p = -1 we shall have 
2 

E[Z + Y] = 0, and Y = -Z with probability one. Thus if both 
bounds are attainable Z must have a symmetric distribution. 

Given general marginal distributions F^(y) and 
Mardia (1970) showed Frechet bounds as 

max [0,F^ (y) tF^ (z) -1] ^F(y,z) _< min [F^^ (y) , F^ ( z) ] (II-A-1) 

From this we can find the range of possible values of p. 

For simplicity we now confine our consideration to distribu- 
tions F(y,z) of positive random variables whose deriva- 
tives F^(y) and F^Cz) are strictly positive for y > 0, z > 0, 
respectively. Suppose also that the variates are scaled to 
have unit variances. Let Gj^(u), G 2 (v) be the inverse func- 
tions of F^(y), F 2 (z), i.e. 

F^[Gj^(w)] = w 

where 0<w<l, i=l, 2. Then the correlation coefficient 
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between Y and Z is given by the equation 



00 



00 



p 



J / y z dF(y,z) - E[Y] E[Z] 
0 0 



(II-A-2a) 



1 1 

/ / G, (u) G (v) dK(u,v) - E[Y] E[Z] 

0 0 -^ 



(II-A-2b) 



where K(u,v) is the joint distribution of the quantities 
U = Fj^(y), V = F 2 (z). Then U,V are jointly distributed on 
the square 0_<U<_1, 0_<V_<1, in such a way that the 
marginal distributions are uniform on the unit intervals. 
From expression (II-A-1) the minimum correlation is attained 
when the probability is concentrated uniformly on the line 
U + V = 1. The minimum value of p then is 



The corresponding F(y,z) will be a singular distribution with 
all the probability concentrated on the line Fj^(y) + = 1. 

In fact Y and Z are an antithetic pair. The maximTom value 
of p is attained when the probability is concentrated uni- 
formly on the line U = V. Then the maximum value of p is 



Pmin 



0 



1 

J G^(u) G2(1 -u) du - E[Z] E[Y] (II-A-3) 



P 



max 



0 



1 

/ G^(u) G 2 (u) du - E[Z] E[Y] 



(II-A-4) 
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with the corresponding probability concentrated on the line 
(y) = * 

By using this result we will get lower and upper bounds 

of the correlation coefficient for uniform, exponential and 

gamma marginal cases. In the uniform marginal distribution 

case, we can see p . = -1 and p =1 from Moran's condition 

min '^max 

ii) , i.e., uniform distribution is symmetric about its mean. 
For the exponential marginal distribution case suppose Y and 
Z have exponential distributions with unit mean and variance. 
Then 



F^(y) = 1 - e~^, = 1 - e"^ 

and the inverse functions are 

(u) = - In (1 - u) , 

^2 (v) = - In ( 1 - v) . 

From equations (II-A-3) and (II-A-4) 

1 

Pmin = G2(1 -x) dx - E[Y] E[Z] 

1 2 
= f In X In (1-x) dx - 1 = 1 - ^ 

0 ^ 

: -0.64493 
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1 

Pmav " I G, (x) G„(x) dx - E[Y] E[Z] 
max „ i i 



1 

= / In X In X dx - 1 = 1 

0 



The p . can be attained when all Y and Z are concentrated 
min 

, — y —2 

on the line e -^ + e =1. Also the p can be attained 

max 

“V “ z 

when all Y and Z are concentrated on the line e = e , 
i.e. , Y = Z . 

For the gamma marginal distribution case, suppose that 
Y and Z have gamma type distributions with probability density 
functions 



^1^^^ " a > 0 , y > 0 

6 > 0 , z ^ 0 . 



Then 



E[Y] = a, E[Z] = 6, VAR[Y] = a VAR[Z] = 0 . 



The will be attained when the probability density is 

concentrated on the line 



1 

r (a) 




0 



a-1 

X 



dx + 



r (B) 



-X 



, 6-1 



dx 



1 . 
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This defines y uniquely as a function of z which can be 
written y = A(z) . The is then, on rescaling the 

covariance 



P • 
min 



{ / u A(u) f, (u) du - aS}/(a3)^'^^ 
0 



When a,B become large, p . tends to -1. And the p 

min max 

will be attained when the probability density is concentrated 
on the line 



1 

r (a) 



-X a-1 J 
j e X dx 

0 



1 

r (B) 



2 1 
( -X a-1 j 
j e X dx 

0 



This defines y as a function of z which can be written 
y = B(z). The p then is, on rescaling the covariance 

III 3.x 



P 

max 



{ J u B(u) f , (u) du - aB}/(aB)^'^^ 
0 



when a = B, p tends to 1. Schmeiser and Ram Lai (1979) 
max 

showed the obtainable correlations between random variables 
Y and Z having gamma marginal distribution with density function 

a . -1 

fj_(x) = [(x/6j_) ^ exp (-x/B j_) / (B j_r (aj_) ] 

for X > 0, > 0, Bj_ > 0, i = 1, 2. 
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The Figures (Il-a) and (Il-b) show the obtainable 
correlation as a function of given = 1 and 5, for 

6 , = 1 . 



I 




a-L 



Figure Il-a. Obtainable correlation as a function of 

for = 1, 3^ = I, by Schmeiser 
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Figure Il-b. Obtainable correlation as a function of a 2 
for = 5, 3^ = 1, by Schmeiser 
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B. GENERAL METHOD REVIEW 



1 . Johnson and Tenenbein's General Method 

A general method of constructing a bivariate dis- 
tribution, whose marginal distribution functions are Fj^(x) 
and t is proposed by Nataf (1962), can be represented 

as follows. 

General Method 

i) Consider any two continuous random variables U 
and V with probability density functin h(u,v). 

ii) Let X' = Hj^(u) and Y' = H 2 (v) , where Hj^(u) and 
H 2 (v) are the cumulative distribution functions 
of U and V, respectively, 

iii) Define 



X = f'^(X') = f"^[H^(u)] 



(II-B-1) 



and 



Y = F2^(Y') = f"^[H2(v)] 



(II-B-2a) 



or 



Y = F2^(1-Y') = F2^[1-H2(v)] (II-B-2b) 
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since X', Y' and 1-Y' are uniformly distributed over the 
interval [0,1], X defined by expression (II-B-1) and Y 
defined either by (II-B-2a) or (II-B-2b) will have a joint 
distribution whose marginal distribution functions are Fj^(x) 
and F 2 (y). The method is probably what one would think of 
first but again involves inverses. Based on this general 
method, Johnson and Tenenbein (1979) developed two proce- 
dures for generating (and also constructing) bivariate dis- 
tributions whose marginal distributions and measures of 
dependence, as given by Kendall's T or the grade correlation 
coefficient p^, are specified. Note that these Kendall's T 
and grade correlation coefficient are not the same as the 
ordinary product moment correlation coefficient p we have 
been considering; these measures are discussed by Kendall 
(1962) and Kruskal (1958) and can be defined as follows. Let 
X and Y be continuous random variables having some joint 
probability density function. Let (X^,Y^^), (X 2 /Y 2 ) and 

(X^rY^) be three independent pairs of observations having 
the same joint density function, then 

T = 2 P[(Xj^-X2) (Y^-Y2) > 0] - 1 

Pg = 6 P [ (Xj^-X2) (Y^-Y^) > 0] - 3 

The first procedure for generating bivariate pairs, called 
the WLC (weighted linear combination) , defines 
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u 



U' 



(II-B-3) 



and 



V 



cU' + (l-c)V 



(II-B-4) 



where 0 ^ c ^ 1, U' and V are independent and identically 
distributed random variables with probability density func- 
tion g{*). Then X and Y are obtained from the general 
method using equation (II-B-1) and (II-B-2a) if the depen- 
dence measure is positive or equations (II-B-l) and (II-B-2b) 
if the dependence measure is negative. 

In order to apply the WLC procedure, we must obtain 
expressions for H^(u) , H 2 (v), a^{u,v) and T(u,v), in terms 
of c and g(*). The values of Hj^(u) and H 2 (v) allow us to 
apply the general method for a given choice of c and g(t). 

The expressions for T (u,v) and (u,v) allow us to specify 
c for a given choice of g{*) in terms of the required value 
of either T or . From equations (II-B-3) and (II-B-4) 



u 



H^(u) 



/ g(t) dt 



(II-B-5) 



— 00 



H2(v) 



/ / g(u') g(v') du'dv' (II-B-6) 
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where 



= {(u',v'); cu'+(l-c)v' ± v} 



and the joint density function of u and v is 



h(u,v) 



1 , . ,v - cu. 

g(u) g(T^) 



(II-B-7) 



00 00 



(u.,v) = 12 I j (u) H 2 (v) h (u, v) dudv - 3 (II-B-8) 



— 00 — 00 



T(u,v) = 4 I G„(^-^t) g (t) dt , (II-B-9) 



where 



92 (t) = / g(w + t) g(w) dw 

— 00 

00 

G 2 (t) = / g 2 (x) dx 

— 00 

Using equations (II-B-5) , (II-B-6) , (II-B-7 ) , (II-B-8) and 

(II-B-9), we have to evaluate Hj^(u), H 2 (v) , h(u,v), T (u,v) , 
and pg(u.v) for all values of c and for specified g(*)* 

As Johnson and Tenenbein mentioned in their report, these 
integrals are quite tedious to perform, so they gave some 
computation results in their report. 
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The second procedure, called TVR (trivariate reduction) 
is discussed by Mardia (1970) . In this case U and V are 
defined as 

U = U' + ez' (II-B-10) 

V = V + BZ' (II-B-11) 

for 0 ^ B < «, and U', V, and Z' are independent and iden- 
tically distributed random variables with probability density 
function g ( * ) . 

Like the WLC procedure, one needs to define H^(u), 

H„(v), h(u,v), p (u,v) and T(u,v) as a function of B and g(*). 

s 

From equations (II-B-10) and (II-B-11) it follows that 

H^(u) = H 2 (v) = / /g(u') g(z') du' dz ' (II-B-12) 



where 



= {(u',z'); u' + Bz' _<u}. 



00 

h(u,v) = / g (u - Bz) g (v - Bz) g(z) dz (II-B-13) 

— 00 



and 
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T (u, v) 



(II-B-14) 



00 

=4 J [G^{Qt)]^ g^{t) dt - 1 

— 00 



where 



92 (t) = / g(w+t) g(w) dw 

— oo 

00 

G2(t) = / 92 (x) 

— 00 

using equations (II-B-12) , (II-B-13) , (II-B-14) and (II-B-8) 

the same computations are needed. 

In this general WLC and TVR methods, there are two 
problems. One is the need for inverse transformations. The 
other is the need to generate coupled random variables to 
create dependence after the inverse transformation is applied. 
Doing this by linear combinations is not necessarily the 
simplest and most felicitous with respect to calculation of 
correlations and range of correlations. 

A significant simplification can be achieved by 
obtaining a "smooth" bivariate uniform pair (U,V) whose 
correlation can be varied between the limit of correlation 
which can be obtained for uniform marginals. These limits 
are (-1,1) since the antithetic pair (U,l-U) has correlation 
- 1 . 
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Note that the problem is essentially one of obtaining 
a bivariate uniform random variable with positive correlation 
because 



corr(U,V) = p 

^u,v 



"Pu,l-V = corr(U,l-V) 
~‘^1-U,V ^ corr(l-U,V) 



Two fairly simple schemes for obtaining bivariate uniform 
distributions are discussed in the next section. 



C. REVIEW OF SOME SPECIFIC METHODS 

The mixture-truncation method is a general algorithm 
for generating bivariate pairs (Y,Z) with given marginal 
distributions and requires only the availability of a generator 
for the marginal distribution and calculation of constants. 
Unless randomization is used, however, it does give a bi- 
variate distribution with a discontinuous joint distribution, 
as in Section (III-D) . This could be a disadvantage. It 
is interesting therefore to compare it to other generations 
which are specific to the given marginal distribution and 
we do this here for the cases for which the mixture-truncation 
method is illustrated in this thesis, namely uniform, exponen- 
tial and gamma. 

There are many schemes available for uniform and exponen- 
tial cases though it should be noted that smooth, simply 
generated pairs with easily calculated correlations have 
only recently been available. In particular we give here two 
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very recent schemes for smooth bivariate uniforms which are 
easily generated and whose correlations can be calculated; 
similarly two schemes for exponential which are essentially 
the uniform schemes after transformation. When one comes 
to gamma marginals one is again in difficulty and we describe 
a recent scheme by Schmeiser and Ram Lai which involves 
however, computation of the inverse gamma distribution 
function and very difficult initialization computations. 

The mixture-truncation method is quite competitive here. In 
fact if one can compute the inverse gamma distribution then 
bivariate gamma's can be computed as [Y = F ^(U), Z = F ^ (V) ] 
where U and V are a bivariate uniform generated by the schemes 
mentioned above; this would be simpler for any marginal dis- 
tributions than some of the general schemes mentioned in the 
previous sections. Comparisons of the generating schemes given 
in this section with the mixture- truncation method are given 
in Section IV-D and Section V-D. 

1 . Lawrance and Lewis * s Bivariate Uniform 

Lawrance and Lewis (1979) show that a simple trans- 
formation of the NEAR(l) (New Exponential Autoregressive) 
process gives a two-parameter family of Markovian random 
variables with uniform marginal distributions. This method 
generates a correlated uniform pair as a multiplicative mix- 
ture of uniform variables, where the parameters a and B take 
on values between zero and one, with the condition that they 
not both be one. Let Y be a uniform [0,1] random variable, 
and define 
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wp a 



e Y 



S 



z 



wp 1-a 



where 



U 



wp 1-e 

1 - (l-a) 6 



U 



(l-a) e 




where U is an independent identically distributed uniform 
(0,1) random variable. The correlation structure is defined 
as follows as a function of a and 3; 



The model can be reduced to a one-parameter model by suitably 
relating a and 3, e.g., a = 3 and all positive values of p 
can be obtained. Consequently Y and 1-Z will have a full 
range of negative correlations. A generating procedure for 
this bivariate pair of uniform random variables is as follows. 

Generating Procedure 

1. (Initialization) . For given correlation, find suitable 
parameter values a and 3 • 



P 



Y,Z 



3 , a 3 

2 + 3 1 + ( 1— ct ) 3 



2. Generate U, , set Y U, ; r = l-a, P ^ 

3. Generate U 2 . 
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4. 



^2 

If U 2 £ r, set U ■«- — and go to 8. 



5. Otherwise U -t- 



U2-r 
1 - r* 



6. If U ^ P, set Z and go to 10. 

7. Otherwise, set Z -e 



8. If U ^ P, set Z ^ and go to 10. 

9. Otherwise, set Z . 

10. Go to 2 and repeat until the desired number of bi- 
variate uniform variables (Y,Z) are obtained. 

The program is listed in the appendix and the scatter plot 
and bivariate histogram of resulting output are in Section 



IV-D. 



2 . Bivariate Uniform By Transformation of Gaver's 
Bivariate Exponential 

In Gaver's bivariate exponential to be discussed in 
Section II-C-4, we can define a bivariate exponential random 
vector, (Y,Z) with correlation p = - P/2 as follows: 



Y 



ln[ 



U^/N (i_p) 



1-P 



(II-C-1) 



and Z is, conditional upon N = n, a gamma random variable 
with shape parameter n and mean n(l-p), where N is a geometric 
random variables with probability density function 
f(x) = p(l-p) , X = 1, 2, ... and U is a uniform random 

variable over [0,1]. From this we can generate a pair of 
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uniform random variables (V,W) by transformation. Since we 
know that to generate an exponential random variable X from 
a uniform random variable U[0,1), we use the following 
transformation: 



X = - ln(l - U) 



(II-C-2) 



The resulting X has the exponential distribution with unit 
mean. From equations (II-C-1) and (II-C-2) , we can generate 
V, a uniform random variable over [0,1]: 



V 



^ (1-p) U 

1 - P U 



1/N 

1/N 



Further also we can generate W, a uniform random variable 
over [0,1] : 



1-P 

W = ( n U. ) 

• 1 1 
1=1 

This follows since we know that Z has, conditionally upon 
N = n, the gamma distribution with shape parameter n and 
mean is n(l-p) . 

In general the resulting V and W will be negatively 
correlated and (V,l-W) will be a positively correlated pair. 
Because the correlation structure will be changed by trans- 
formation, the resulting (V,W) need not have the same 
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correlation with (Y,Z). The correlation coefficient of 
(V,W) will be computed as a function of P as follows. 



p = 12 E[VW] - 3 



where 



E[VW] = Ej^ [E [VW In = n] ] 




ilzPlP 

1 - P U 



1/n 

T7K 



du 



[ du]''] 

0 



Unfortunately this computation of p as a function of P is 
difficult. As an example we used here the same function for 
p as holds in the exponential case. 



Generating Procedure 
1. (Initialization) 

i) Compute P for given correlation p 
ii) Choose N from a geometric distribution with 
parameter 1-P 



2 . 



Generate a uniform [0,1] random variable and define 



V = 



(1-P) U. 



1 - P U. 



1/n 

TAT 



3. Generate N = n uniform [0,1] random variables 
i = 1, ..., n, and define 
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w = ( n u. ) 

i=l ^ 

4. Deliver (V,W) and go to 2 until a sufficient number 
of bivariate pairs have been generated. 

The program is listed in the appendix and the scatter plot 
and bivariate histogram of resulting random vectors are in 
Section IV-D. 

3 . Marshall and Olkin's Positively Correlated 
Bivariate Exponential 

Suppose X is the age, or length of service of the 
first device at the time of death, governed by two indepen- 
dent Poisson processes Z^{t) , with parameters 

respectively and Y the same of the second device, 
governed by two independent Poisson processes Z^^t.) , 
with parameters X^^ 2 ' respectively. Further suppose the 

second device is placed in operating after a time 6 later 
than the first device. Then the joint distribution of (X,Y) 
is defined as 

P [X > X, Y > y] s F (x,y) 

• exp{-X^x-X2y - 3 ^ 6^0 

- exp{-X j^x-X 2 y~X^ 2 ’^^^ [y (Y' ^ 6^0 

(II-C-3) 

If 6 = 0, then equation (II-C-3) reduces to 
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P [X > x,y > y] 



F (x,y) 



= exp{-Aj^x - max(x,y)} 

Note that since the marginal distributions do not depend on 

6, 



E[X] = 






VAR[X] = 



(X^+X^2^ 



2 ' 



E[Y] = 



^2 ■^^12' 



VAR[Y] = 



(X 2 + Xj^2^ 



We also have, for 6 > 0 



E[XY] = 



^12 ® 






^^l''’^12^^^2''’^12^ ^^l'*'^12^^^2''’^12^ 



and the correlation 



corr(X,Y) = p 



x,y 



Xj^2 “ ( A j^+X ^^ 2 ) ^ 



(II-C-4a) 



where 



X X ^ X 2 X ^ 2 ' 
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If 6 = 0 , then the equation (II-C- 4 a) reduces to 



P 



xy 




(II-C- 4 b) 



This characterization can be represented in terms of 
independent random variables. For 6 ^ 0 , if there exist 
independent exponential random variables and 

with respective parameters \ 2 > ^12' ^12' 



X = min(U^,U2) 



minCU^/U^-S) 


if 


U3 > 6 


. min(U2,U^) 


if 


U3 < 5 



It can be verified formally from the relation 



P[X>x,Y>y] = P > X, U2 > y ] {P [U^ > x,!;^ > y +<5 1 > iS ] 
X P[U^ > 6]} + P[U^ > y]P[U^ > xlu^ < 6]P[U2 > 6] } 

In the 6=0 case, the representation yields 



X = min(U^,U2) 
Y = min(U2,U2) 
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This bivariate distribution has a line-discontinuity, if 
E[X] = E[Y], along the line x = y, i.e., in the case where 
X = y = U^. Thus this bivariate exponential is not as smooth 
as others. In particular the NEAR(l) process of Lawrance 
and Lewis generates bivariate exponentials with a continuous 
density function. But if we consider the case 6^0, the 
line-discontinuity will be removed. In either case, the 
resulting (X,Y) pairs always have positive correlation as 
shown in equations (II-C-4a) and (II-C-4b) . If we use these 
methods to generate positively correlated random vectors 
(X,Y) for exponential marginal distribution with unit mean 
and given p, we have to compute X^, X 2 and X ^^2 given 

correlations to generate independent exponential random 
variables. For simplicity we will consider the method with 
6=0. From the relationship X = Xj^ + X 2 + '^■^ 2 ' 



E[X] 



'1^ 



^12 



1 



E[Y] 



X2 + 



12 



P 




1 



we can compute X^^, X 2 and X ^^2 functions of p. 

X j^2 “ (1 + p ) 



X 



1 




1 



X 



12 • 
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Generating Procedure 



1. (Initialization) . Compute and for given 

0 ^ p ^ 1 • 

2. Generate three independent exponential random variables 

^1' ^2 ^3 unit mean. 

3. Rescale the independent random variables and set as 



Ui ^ 



Uj * E/Xj 






4. Define X and Y as 



X = min[U^,U2l 
Y = min[U2,U2] 

5. Deliver (X,Y) and go to 2 until a sufficient number 
of bivariate pairs have been generated. 

The program is listed in the appendix and the scatter plot 
and bivariate histogram of resulting random vectors are 
shown in Section V-D. 
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4 . Gaver ' s Negatively Correlated Bivariate Exponential 
This negatively correlated bivariate exponential 
is generated from the following considerations. Suppose a 
particular system has N defective elements; here we assume 
N has geometric distribution with probability density 
function 



The time to failure of each defective element is T. with 



density function F(t). Further suppose repair is carried 
out after all defective elements are discovered. If we 



provided N = n, and define T as minimum of T^^, then the 
joint distribution function of T and R is obtainable from 



f (x) = p(l - p) 



x-1 



X = 1, 2 



1 



define R as repairing time which is distributed as G^(x) 



P [T > t,R <_ x] = [ (1 - F(t) ) 

n=l 




where 



(1 - p) P 



n-1 



0<p<l, n=l, 2 



P 



n 



Then 



E [R|T > t] 



E[R] 



1 - P[1 - F(t) ] 
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where E[R] is the expected value of an element repair time. 
Moreover, 



E[r|t = t] 




E[R] [1 + 3^(1 - (|.(t) ) ] 



where 



4>(t) 



P [T < t] 



F(t) 



1 - P[1 - F(t) ] 



From this regression, we know that if T is long, R is short 
and vice versa; the negative relationship is clearly present. 
Further we can calculate Cov[R, T) as 



where T(l) has the distribution of the smallest of a sample 
of two from <J)(t). 



It will not be shown that one can select F(t) = l-F(t) 



in such a way as to induce exponentially distributed time 
to system failure, T. Since we know that exponential G, 
geometrically compounded, yields exponential R, the outcome 
will be a (T^R) pair with exponential marginals and negative 
correlation. 



Cov[R,T) 



P E[R] E[T] (1 - ^ 



< 



0 
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If 



F(t) 



1 - F(t) 



satisfies 



I [F(t) (1-p) 



e 



-xt 



n=l 



then T is exponential with mean 1/X. Solution for F 
yields , 



which is a logistic distribution, left truncated at t = 0. 
If X = 1, then T is exponential with unit mean, while 
E ] = 1/2 , so 

corr(R,T) = - £ E[R] E[T] . 

If G is chosen so as to make R exponential then it may be 
shown that 



Consequently, a greatest lower bound for the correlation in 



F(t) 



1 



t > 0 



p + (1-p) e^^ 



corr (R,T) 




this model is - j. 
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Generating Procedure 

1. If 0 < p < -0.5, Set P -e -2p 

2. Generate a geometric random variable N with parameter 1-P 

3. Generate a uniform (0,1) random variable U and define 



4. Generate a gamma random variable, Z, with shape 
parameter n, and mean is n(l-p) 

5. Deliver (Y,Z) and go to 2 until a sufficient number 
of pairs are obtained. 

For obtaining negatively correlated bivariate exponential 
random vectors, this algorithm is very simple and one of the 
few available. Moreover the correlation is known explicitly. 
The program is listed in the appendix and scatter plot and 
bivariate histogram of resulting random vectors are shown 
in Section V-D. 

5. Arnold's Vibariate Gamma Generator 



to generate bivariate gamma random vectors having positive 
correlation with 



Y as 



Y 





Arnold (1967) developed a Trivariate reduction method 
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where and are the given shape parameters for the marginal 
distributions. Letting "gainma (a,g)" denote the gamma dis- 
tribution with shape parameter a and scale parameter g , 
so that 



- FrTST ' 

where r{a) is the gamma function and 

E[x] = ag , VAR[x] = ag^. 



the trivariate reduction algorithm proceeds as follows for 
given shape parameter values 0.2 and p to generate unit 
scale gamma variate i.e., g = 1 . 

Generating Procedure 



1. 


Generate 


^1 ' 


gamma 


- p ( 


a/2 

03^02) /i) 


2 . 


Generate 


G2 , 


gamma 


{ 0-2 - P ( 


a/2 

03^02) /i) 


3 . 


Generate 


^3 ' 


gamma 


(p (03^02 




4 . 


Define Y 


and 


Z 







Y 



Gi + G3 



Z 



G2 + G3 
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These Y and Z are unit gainina variates with 6 = 1, but Y and 
Z can be multiplied by 6 j^ and 62 » respectively, to obtain 
any desired scaling. In the algorithm, is a common com- 
ponent of both Y and Z, thus inducing positive correlation. 

When a positive, but not extreme, correlation is 
needed, trivariate reduction is an excellent algorithm, 
using univariate gamma generators. If the marginal distri- 
butions have the same shape parameter value a, then the 
correlation limitation is 0 ^ p < 1. But in the other case 
the upper limit is smaller than 1 . 

Note that this algorithm exploits the infinite 
divisibility of the marginal gamma variables and is applicable 
to any infinitely divisible marginal. It is commonly used, 
for instance, to get bivariate Poisson random variables. 

6 . Schmeiser's Bivariate Gamma Generator 

Schmeiser (1979) developed a family of algorithms, 
any member of which can generate bivariate gamma random 
vectors having any shape parameters a^, and allowable 
correlation p. In this section we will discuss his general 
procedure . 

If Z, Wj^ and W 2 are independent gamma random varia- 
bles with shape parameters r, 6 ^^, and respectively, and 

U is an independent uniform [0,1] random variable, and either 
V = U or V = 1-U; then 

X, = f“^(U) + Z + Wt (II-C-5) 

In 1 
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^2 " + Z + «2 

where F ^ and F ^ are the inverse distribution functions 
n m 

of gamma (n,l) and gamma (m,l), respectively, distribution 
function. The resulting X^, X 2 are gamma random variables, 
with shape parameters 



= n + r + 5^, 

(II-C-6) 

“2 ~ m + r + 



respectively. This result follows immediately from the 
reproducibility property of the gamma distribution and from 
noting that F ^(u) and F~^(l-u) are each random variables 
having CDF F. This scheme generalizes the trivariate 
scheme but brings in the inverse CDF. 

The correlation coefficient is defined as 



E[f"^(u) 

n 



.-1 

m 



(v) ] - nm + r 



Xi,X2 






1/2 



(II-C-7) 



The remaining problem is to select values of the five 
parameters n, m, r, 6^ and <S 2 to obtain the desired marginal 
distributions and correlation. The following conditions 
must be satisfied. 



n + r + 
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m + r + 6^ = 



E[F ^(u) F ^(v)] - (nm) + r 
n m 






1/2 



= P 



(II-C-8) 



1 ^ 2 ~ ^ 



Since we are using five variables to satisfy three equality 
conditions, finding a set of parameter values corresponds 
to finding a feasible solution, rather than an optimal solu- 
tion, to a nonlinear programming problem. 

An efficient solution procedure for determining 
parameter values is important, since substantial computation 
is required to determine whether or not conditions (II-C-8) 
are satisfied for given parameter values. Most of the 
computation is involved in calculating 



p 



E[F"^(u) f"^(v)] 
n m 



(a^a2) 



1/2 



nm - r 
' / 



since the expected value must be calculated numerically 
using any one of the following three integrals: 



/ f“^(u) f“^(u) du (II-C-9a) 

0 ^ n m 



0 




(F^(x)) 



x^ exp(-x) dx/r(n) 



(II-C-9b) 
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(II-C-9C) 



/ F^^(F^(-ln y) ) (-In y) ” dy/r (n) 

if p > 0. If p < 0, then replace F^^(u) with F^^(l-u) in 
equation (II-C-9a) , replace with 1 - F^(x) in equation 

(II-C-9b) , and replace F^(-ln y) with l-F^(-ln y) in equation 
(II-C-9C) . 

Schmeiser seemed to have best results in terms of 
a subjective trade off between speed and accuracy, using a 
24 point Gaussian method for the integration. He selected 
the parameter values from the feasible values satisfying 
conditions (II-C-8) by making the curves of regression 
E[Xj^|X 2 ] and E[X 2 |Xj^] behave as desired. 



Generating Procedure 



1. 


Generate X^ ~ gamma 


(n,l) 


2. 


U ^ F^(X, ) 
n 1 






If p < 0, U^l-u 




3. 


Generate Z - gamma 


(r,l) 


4. 


Generate - gamma 


(6j_,l) 


5. 


Generate W 2 ~ gamma 


1 — 1 
<N 


6. 


Xj^ Xj^ + Z + 





X 2 ^ f"^(U) + Z + W 2 
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Based on these procedures he developed a family of algorithms 
which can provide variates having any theoretically possible 
correlation p. Anyway, for gamma marginal distributions, 
not all correlations are consistent with particular shape 
parameter values. Schmeiser shows the obtainable correlations 
as a function of given = 1 and 5 as shown in Section 

Il-A. Note that only when is it possible to obtain 

p = 1- Likewise p = -1 Is not possible except in the limit 
as and tend to infinity. The maximum, and m.inimum 
possible correlations, given in Moran (1969), occur when 



distribution function and inverse CDF, respectively, of the 
gamma distribution with shape parameter a. 



X 



2 




and 





respectively, where F (X) and F ^(U) are the cumulative 

a a 
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III. GENERAL MIXTURE-TRUNCATION METHOD 



Denote by (Y,Z) the bivariate random pair, where each 
has identical marginal continuous distribution F(x), and 
denote a general random variable from this distribution by 
X. The argxment is not specific to continuous random varia- 
bles; this aspect comes in only in the computation of the 
correlations and can be developed in a parallel fashion for 
discrete marginal distributions. Let 



P = 



1 - 



I-CL 2 ^2 



with stationary vector 



1 - a. 



1 - a. 



IT = 7T P = ( 



1 — 02^ + 1 — 02^ 1 — 0j^+l-02 






and let x, be an X truncated to the left of a fixed point x , 
i o 

X be an X truncated to the right of x , so that 
^ o 



Fj^^(x) = P[X^<x] 



F(x) 

F(x^) 



if X < x_ 

— o 



if X > X 
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0 



if X < X 



o 



F (x) 
x„ 



2 



P[X2 ix] 



F(x) - F(x^) 



if X > X 



1 - F(x^) 



o 



In addition we set tt^ = F(x^), 



1 - 7 T^ and choose 



Y and Z as follows. 

i) Choose Y from with probability and then choose 
Z from with probability or from X2 with 

probability 1 - a^. 

ii) Choose Y from X wi-l-h rvrnhah-i l i -f-v TT _ wh TT_+Tr = 3.^ 




If we choose (Y,Z) as in the above procedure, then we can 
make the following two theorems. 

A. MARGINAL DISTRIBUTIONS 
Theorem 1 . 

The marginal distribution of (Y,Z) becomes F(x) for both 
Y and Z. 

Proof 

1 . Marginal distribution of Y 

By definition Y is the mixture of and X2 with proba- 
bility ir^, Tt^ respectively. That is 




and then choose Z from X^ with probability 1-02/ or 
from X2 with probability 02* 




TT^ (x) + 1^2 (x) 
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F^(x) = 



"l + "2 • 0 



F(x) - F(x^) 

""l * ^ ^"2 1 - F(x^) 



if X < X 

— o 



if X > X 



But since we define = F(x^); tt 2 = 1 - = l-F(x^), we 

have 



Fy(x) = 



F(x ) 

.1-^ = 



F(x) - IT. 



TT, + 7T 



12 TT. 



if X < X 

— o 



= F(x) if X > X, 



= F(x) in all cases. 



So Y has the marginal distribution F(x) . 

2. Marginal distribution of Z 
If Y is from then 



, (x) = F (x) + (l-ai)F (x) 

‘1 ^ ^1 ^2 



F (x) a. / 1 T • n 
‘^1 F(x^) ^ 



if X < X 

— o 




1 + (1 - a^} 



F(x) - F(x^) 
1 - F(x^) 



if X > X 

o 
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If Y is from X^, then 



F_ (X) = (1-a ) F„ (X) + a, F^ (x) 

^2 Z Z 



“2> fiir + “2 ■ “ 



F(x) - F(x^) 
(l-“2> • 1 + “2 ' T - F U") ~ 



if 



if 



So, 



F 2 (x) = TT^ F^ (x) + TT^ F^ (x) 



F(x) , , F(x) 

■'T-i a, -- + (1 - a«) 



1^1 F(x ) "2'-" "2' Fix ) 

o o 



Fix) - F(Xq) 
7r^(a^ + (1 - a^) 1 - F(x^) 



if 



if 



F(x) - F(x^) 
+ n^iil-a^) + a2 i i f ' CxT 



and we defined and tt 2 as follows. 

^ ~ °^2 

^1 1 - + 1 - 02 

^-“1 

^2 1 - + 1 - 



X < X 

~ o 



X > X 

o 



X < X 

— o 



X > X 

o 
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From this, we know 



1 1 - a_ = -^ (1 - a,) 

2 tt2 1 

1 

If we use this relationship, then ~ ^ both 

cases. So Z also has the marginal distribution F(x). The 
I result is a consequence of the fact that it is defined to be 

I 

: the stationary vector associated with P. 

I B. THE PRODUCT-MOMENT CORRELATION 
Theorem 2 . 

The correlation coefficient between Y and Z becomes 

P = 6 M , 



where 



-1 £ 6 = - (1 - a 2 > £ 1 , 



and 



M 





TT, Tr_/a 
1 2' X 



/ 



where 



X 



X d 



F (x) 

P(Xo> 
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oo 



F(x) - F(x ) 

J X d^i =T7 — 

^ 1 - F (x ) 

X o 

o 



J F(x^) 



- y . 



2 ^ F(x) - F(x^) 



I X d ■ 

x_ 



F(x ) " ^2 



X 



/ x^ dF(x) - E[X] 



Proof 



2 2 2 
= 0^ °2 ^2 ” ’^ 2 ^ ^1 ^2 



_ cov [Y y Z ] 



'YZ 






E [YZ] - E [Y] E [Z] 
Oy 



Now 
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E [YZ] 



where 

Further, 

E[Y] = 

by Theorem 



E^[E[YZ|Y,Z e S] ] 

W 

E[YZ|Y e X^,Z e Xj^]P[Y e X^,Z e X^^] 

+ E[YZ1y e X^,Z £ X2]P[YeXj^,Z e X 2 ] 

+ E[YZ[Y£X 2 ,Z £ Xj_]P[Y £ X 2 ,Z £ X^^ ] 

+ E[YZ|Y £ Y.^,Z £ Xj^]P[Y £ X 2 ,Z £ X 2 ] 

E[X^] E[Xj^] Ti^a^ + E[Xj^] E[X 2 ] Tr^^d-a^^) 

+ E[X 2 ] E[Xj^] v^n-a^) + E[X^] E[X^] 

2 

Ml cii TTi + y2(l-ai) tTi 

+ y^ y2(l-a2)^2 ^ ^2^ “2^2 



yi = E[Xi], y2 = E[X2l 



Eg[E[Y|Y £ S]] 

E[Y|Y £ Xj_]P[Y £ X^] + E[Y|Y £ X 2 ]P[Y £ X 2 ] 
E[Xi] Ti^ + E[X2] Ti2 
^1 ^2 ^2 ~ E[Z] 
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Also by Theorem 1, 



2 

ay 



If we put together these formulae into 

E[YZ] - E[Y]E[Z] 

^Y,Z Oy 



E[Y^] - (E[Y])^ 



2 _ 2 



we get 



'Y,Z 



^^1 ~ ^2^ tTiTT 2 (oti - (1 - 02 ) ) 



X 



and let 



6 = - (1 - 0.2) 



M 




2 - 2 
^1^2/‘^X 



then 



p = 6 M . 

C. GENERAL ALGORITHM 

We give here three algorithms for implementing the 
bivariate mixture- truncation method, which we call the FXO 
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I 



method, the UXO method, and the TXO method. All of these 
' methods are exactly the same except in how the algorithm 
chooses x^, the truncation point, from the x^ range [x^,x^] . 
The first procedure, called the FXO, choose x as a fixed 

I 

point of x^ and x^ and uses the same x^ during the entire 

I routine. The second procedure, called the UXO, chooses x 
! ® 

uniformly from ^nd repeats this step in every call 

to the algorithm. The third procedure, called the TXO, is 
: the same as the UXO procedure except in that it uses a tri- 
1 angular distribution instead of uniform. It is necessary to 

i 

I fix these choices of x^ becuase in general there is more than 
j one x^ which will give a bivariate pair (Y,Z) with the given 
I marginal distribution and given correlation. The first 
procedure, FXO, is defective in terms of their discontinuity 

! 

of distribution while the second and the third, UXO and TXO, 
are satisfactory in this respect. The choice of the midpoint 
of the interval for x^ in FXO is based on experience 

presented later for various marginal distributions. Note 
that the algorithm described here is inefficient in that it 
generates the truncated variables X^^ and X 2 by comparing ran- 
dom variables X to x^ until one which is respectively greater 

than or less than x is found. More efficient methods can 

o 

be found in special cases such as the exponential, but the 
present algorithm requires only a generation of univariate 
random variables X without regard to the method used to do 
this. Of course initialization is required and this is 
specific to each marginal distribution. 
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General Mixture-Truncation Method 



1. (Initialization) 

i) For given marginal distribution F(x) and correla- 
tion coefficient p find x^ ranges [x^,x^] 



2. Define truncation point x 



* FXO method 



i) 



2 



* UXO method 

i) Generate a uniform [0,1] random variable 

ii) X = X. + (x - X, ) * 






u i. 



* TXO method 

i) Generate two uniform [0,1] random variables 

Vi, V^. 

^ + X^ -t- X2 

where 



X = ^(x + X ) 

m 2 2 , u 



* ''l 



X- = (x - X„) * V- 

2 urn 2 



3. Compute parameters value, ir^, 02 



4 . Choose type for Y 

i) Generate a uniform [0,1] random variable U 
ii) If U £ TTj^, go to 9 
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5. Y is an 

i) Generate a random variable X from F(x) 

ii) If X > x^, set y ■<- X and go to 6 

iii) Otherwise return to 5.i) 

6. Choose type for Z 

i) Set U ^ ( (U - TT^)/(1 - TT^) ) 

ii) If U _< 1-02/ go to 8 

7. Z is an X 2 

i) Generate a random variable X from F(x) 

ii) If X > X , set Z X and go to 11 
o 

iii) Otherwise return to 7.i) 

8. Z is on Xj^ 

i) Generate a random variable X from F(x) 

ii) If X _< x^, set Z X and go to 11 

iii) Otherwise return to 8.i) 

9. y is an 

i) Generate a random variable X from F(x) 
ii) If X _< x^, set Y X and go to 10 

iii) Otherwise return to 9.i) 

10. Choose type for Z 
i) Set U 
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ii) If U _< a^, go to 8 
iii) Otherwise go to 7 

11. Deliver (Y,Z) and go to 4 for the FXO method, or go 
to 2 for the UXO and TXO method until a sufficient 
number of random vectors are obtained. 

D. BIVARIATE DISTRIBUTION FUNCTIONS 

From Theorem 1, in Section III-A, we know that if Y 
is from X^, then 






“l ^X 



F(z) 



if z < X 




o 




F(z) - F(x^) 



if z > X 



o 



and if Y is from X^, then 



F^{z\y) = ( 1 - 02 ) “2 




if z < X 



o 




F(z) - F(x^) 



if z > X 



o 
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f 



I 

t 

iBy using these we can define the bivariate distribution 
function as follows. 



F(y,z) = P [Y _< y,Z _< z] 



I P [Z <_ z I Y 



u] dP [Y < u] 



where 



P[Y < u] = F(u) 



P [Z _< z I Y = u] 



F(z) 

1 F(x^) 



ct^ + (1 - a^) 



F(z) - F(x^) 

-nr^Fu^ 



(1 - olJ 



F(Z) 

2' F(x^) 



F(z) - F(Xq) 

(1 - a^) + 1 _ F(x ) 



So, if we put these together, integrating with 



if u < X , z < X 

— o' — o 

if u < X , z > x^ 

— o o 

if u > X , z < X 
o — o 

if u > X , z > X 
o o 

respect to 



dP [Y ^ u] , we get the final result: 
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P[Y <_y,Z_< z] 



For example 



- ^F(z)F(y) ify^x^, z^x^ (III-D-1) 

(1-aJ 

■ ^” 1 + [F(z)-F(x^) ]}F(y) if y£X^, z > x^ (III-D-2) 

(1-a, ) 

- {a, + ^[F(y)-F(x )]}F(z) if y>x^, z < x^ (III-D-3) 

X ^2 o o — o 

a^TT^ + (l-aj^) [F(z)-F(x^)] (III-D-4) 

“2 

+{ ( 1 - 02 ) +:jp[F(z)-F(x^) ] } [F(y)-F(x^) ] ify>x^, ^ ^ 
the expression (III-D-4) is obtained as 



F(y,z) = / P[Z <_ z|y = u] dF(u) 



^o' 




J P[Z ^ z|y = u ^ X ] dF(u) 

— 00 



+ I P[Z^z[y = u> x^] dF(u) 

^o 

(1-a,) 

= {a, + — [F(z) - F(x )]} F(x ) 

1 ^2 o o 

+ (l-a_) + — [F(z) - F(x^)]}[F(y) - F(x^)] 

2 T72 o o 

It is easily seen from (III-D-2) that when z 
P [Y _< y,Z 00 ] = F(y) ; from (III-D-3) that as y -► »/ 
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P [Y £ _< z] = F(z) and from (III-D- 4 ) that as y 

jP[Y ^ “>/Z ^ z] = F(z) and that z P [Y ^ y,Z ^ «] = F(y) . 

In particular from (III-D- 4 ) we have that, as y >=°, 



F(y,z) a^F{x^) + (1 - tt2) [F (z) - F(x^) ] + {l-a2)tr2 

+ oi^ [F(z) - F(x^) ] 



= aj^F(x^) + (l-a 2 )TT 2 tF(z) - FCxq)] 



= F(z) + ( 1 - 



“ 2'”2 



- ( 1 - 






F(z) 



[where at the last step we used the facts that F(x^) = 
and TTj^{l-aj^) = 712(1-012). If F(x) is absolutely continuous 
■with probability density function f(x) then the joint p.d.f. 
for the bivariate pair (Y,Z) is 



f (Y/Z) 



a, (1 - a, ) 

1-a 

- ^ f(y) f(z) 

l-ai 

- -= — - f(y) f(z) 

^2 

“2 

f(y) f(z) 

^"2 



if y < X , z < X 
^ — o — o 



if y < x^, 



if y > x^, z < x^ 
if y > x^, z > x^ 



(III-D- 5 ) 

(III-D- 6 ) 

(III-D- 7 ) 

(III-D- 8 ) 
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Note that there is a discontinuity in the density function 
as one crosses the boundaries of the four quadrants defined 
by the lines y = x^; z = x^. The density is the same in the 
first and third quadrants. The multipliers of f(y)f(z)/TT2 
are the same in all four quadrants iff there is independence. 
This occurs when ( 1 -a^) = «2 so that = 02 and 
f(y/Z) = f(y)f(z) for the whole range of y and z. 
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IV. BIVARIATE UNIFORM GENERATOR 



The uniform random variable is a continuous random 
variable with probability density function which is constant 
over the interval (a,b) and zero otherwise; the density 
function 



f (X) = 



b-a' 



if a < X < b , 



' 0; otherwise ; 



mean = E[X] = 



b+a 
2 ' 



variance = VAR[X] = 



(b - a) 
12 



We note that if U has a uniform distribution over [0,1] 
then X = a + (b-a)U has a uniform distribution over [a,b] . 

So we will only be concerned with algorithms for generating 
uniform [0,1] distributions. 

Two algorithms for generating uniform bivariate pairs 
are given in Section II, Gaver's and Lawrance and Lewis's, 
and since they have relatively smooth distributed functions and 
are simple to generate, they are probably preferable to the 
uniform bivariate random variable generated by the mixture- 
truncation method unless x^ is taken to be random (of course 
the correlation for the Gaver bivariate uniform is difficult 
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to compute) . However the development in this section of 
the uniform mixture-truncation method does help to fix ideas 
on the problem of initialization and determination of the 
possible range of correlations attainable in other cases. 

To use the mixture-truncation method for generating bivariate 
random vector (Y,Z) whose marginal distribution is uniform 
over [0,1] identically and has given correlation p, we should 
use iiniform [0,1] distribution functions as F(x). Then X^, 
which is X constrained to be on the left side of the trunca- 
tion point x^, becomes uniform over [0,x^] and X 2 / which is 
constrained to be on the right side of x becomes uniform 
over [x^, 1] . 

A. DETERMINATION OF PARAMETERS IN THE UNIFORM MIXTURE- 
TRUNCATION METHOD 

Because X, is uniform (0,x ] and x« is uniform [x ,1] , 

1 ' o 2 " o •' 

E[X^] = = 

VAR[X^] = 

E[X2] = ^2 ^ 

VAR[X2] = 



X 

o 

2 



X 



12 



1 + X 



(1-x^) 

12 



and 



1 - a. 



= = 1-C, tl-c. 



= ==0 ' 



(IV-A-la) 



66 



1 



- a . 



TT- = 1 - F(X ) 

2 . O 



1 “ + 1 - a2 



1 - X (IV-A-lb) 
o 



If we use these formulas in Theorem 2 of Section III, we get 



M = 3 X (1 - X ) 

o o 




p = 3 X (a, - X ) (IV-A-2) 

o 1 o 



» 

iThen from (IV-A-2) , (IV-A-la) and (IV-A-lb) 

I 



+ X (IV-A-3a) 

1 3x o 



a- = 1 - -i(l - ) 

X 2 ^ X 

= 1 ^ 2 (IV-A-3b) 

1 - X 3 ( 1 - X ) 

o o 



Furthermore, and 02 are probabilities so they have to 
satisfy the following conditions: 



0 ^ 1 , 

® 1 “2 — ^ . 
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By solving these two inequality equations, we can get the 
feasible range of for given p as follows. 

x„ = Lower bound of x 
l, o 

r 1/2 - /1/4 -(1/3) p if p > 0 

= (IV-A-4) 

^ /- p/3 if p < 0 

X = Upper bound of x 
u o 

r 1/2 + /1/4 - a/3) p if p > 0 

= (IV-A-5) 

L 1 - /-p/3 if p < 0 

After finding x. and x , choose x within the range (x.,x ] 

by way of either the FXO, UXO, or TXO method. And by formulas 

(IV-A-la) , (IV-A-lb) , (IV-A-3a) and (IV-A-3b) , we can get 

tt 2 / and 02 * Also we can expect some limitations in 

p because of x . From 
o 

M = 3 X (1 - X ) 

o o 

-1 ^8 = - (1 - 02 ) ^ 1 



p = 6 M 
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We can compute that 



p 



max 



3/4 if 6 = 1, <_ = 1/2 j< x^ 



P • 
min 



-3/4 if 6 = -1, X < X = 1/2 < X 

£ — o — u 



In another way, note that if one chooses x^ = 1/2 and 

= 1, from (IV-A-2) we attain the highest allowable 
correlation p = 3/4, and when = 0, p = -3/4, the lowest 
attainable correlation. 

It turns out that for this choice of x = 1/2 (in the FXO 

o 

method) , the bivariate uniform distribution is also about 
the smoothest, as can be seen in a later section (IV-D), 
and from the bivariate distribution. 

B. GENERATING PROCEDURE 

We developed here all of the three procedures , the FXO 
method, the UXO method, the TXO method for generating 
bivariate random vectors whose marginal distributions are 
uniform over [0,1] and correlation coefficient is p . In 
this algorithm we used direct generating procedures for X^ 
and X 2 instead of comparing random variables X to x^ until 
one which is respectively greater than or less than x^ is 
found. 
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Uniform Mixture-Truncation Method 



1. (Initialization) 

i) For given -3/4 < p < 3/4, find x and x 

r 1/2 - /1/4 - 1/3 p if p > 0 

v^- p/3 if p < 0 

_r 1/2 + /1/4 - 1/3 p if p > 0 

X “ 

u 

1 - /^^73 if p < 0 



2. Define truncation point x^. 



* FXO method 



i) X, 



( x , + x ^)/2 



* UXO method 

i) Generate a uniform [0,1] random variable 

ii) X = X + (x - X ) * U, 
o 2 , ' u z 1 



TXO 


method 


i) 


Generate 




Vi, V2 


ii) 


X = X 

o 




where 




^1 






m Z. 
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x_ = (x - x„) * V- 

2 u m 2 






3. Compute and a 2 as 



"l ' “ ==0 



"2 ' ^ • "l 



^ X + X 

3 o o 



1 - -^(1 - a ) 
tt2 1 



4. Choose type for Y 

i) Generate a uniform [0,1] random variable U 



ii) If U ^ TT^ go to 9 



5. Y is on X 2 

i) Generate a uniform [0,1] random variable 



ii) Y -t- X + (1.0 - X ) * W, 
o o i 



6. Choose type for Z 

i) Set U ^ ( (U - TT^)/(1 - TT^) ) 



ii) If U < 1 - a- , go to 8 
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7. 



Z is an X 2 

i) Generate a uniform [0,1] random variable 

ii) Z X + (1.0 - X ) * W_ 
o o 2 

iii) Go to 11 

8. Z is an 

i) Generate a uniform [0,1] random variable W 2 

ii) Z x^ * V?2 

iii) Go to 11 

9 . Y is an Xj^ 

i) Generate a uniform [0,1] random variable 

ii) Y X * W, 
o 1 

10. Choose type for Z 
i) Set U U/tt^ 

ii) If U <_ go to 8 

iii) Otherwise go to 7 

11. Deliver (Y,Z) and go to 4, for the FXO method, or go 
to 2 for the UXO method and the TXO method until a 
sufficient number of random vectors are obtained. 

The programs are listed in the appendix and the scatter 
)lot and bivariate histogram of resulting random vectors are 
shown in Section IV-D. 
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c. REGRESSION OF Z ON Y FOR GIVEN p 



Schmeiser (19 79) in developing a bivariate gaitiina method 

fixed the free parameter by specifying the regression of Z 

and Y and vice-versa. We do this now for the uniform case 

with fixed x and uniformly distributed x . First in fixed 

X case we have, for y < x 
o — o 



These are the step function and clearly not linear in y. 

If X has uniform distribution over [x.,x ], we have the 
w ,x» u 

following results. These results are dependent on Y values. 
If y _< x^, then we have 



E [Z I Y = y /Y IXq] 



E[Xj^] + (1 - a^) E[X2] 




and for y > x 

o 



E[Z|Y = y,y > x^] = (1 - a 2 ) E[X^] + E[X 2 l 




E[Z|Y = y] = / E [Z I Y = y,X = x^,y _< x^l f (x^) dx^ 




— p) dx 

6x X - x„ o 

O U Z 



X 



£ 
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E[z |Y = y] 



(i X - y X - § In — ) 



If x« < y < X , then we have 

X/ — — \Ji 



E[Z|Y = y] = / E[Z|Y = y, X=X^, y>x^]f(x^) dx^ 






X 



u 



+ J E[z|Y = y, X=x , yix ) f(x ) dx 

y 



I 



2 ' 6(T^nrT> <*==0 

X.^ o 



X 



^ 1 

+ J (J 

y 



f (X ) dx 
6x o o 



1 ,1 1 p , 1-y , 

x-x2 u 2 £ 6 y 1-x 

n £ £ 



If y > X / then we 



X 

u 

E[Z|Y = y] = j E[Z|Y = y, X=x^, y > 

X^ 



/ 



■ ^ r -^-pz p) f(x ) dx 

^ 2 6(1- x^) ' ' o' o 



1 ,1 1 p , ^ 

X - ^ X, - ^ In n „ ) 



x-x 2 u 2 Jl 6 1-x 

u £ £ 
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From the results we can see that if we use uniform distri- 

t 

I 

i bution for we can get smoother regression functions if 
I £ y 1. although we still have step functions for the 
jy < x„ and Y > X cases. 

[ J — U 

D. SIMULATION RESULTS 

We will show here the scatter plots and the bivariate 
histograms of the mixture-truncation bivariate uniform (0,1) 

I random vectors with correlation p = 0.3 and p = -0.3 by the 
1 FXO, UXO and TXO methods in Figures (IV-a), (IV-b) and 
I (IV-c) , respectively. 

In addition, we will show the scatter plot and bivariate 
histogram of sample size 2000 from Caver's transformation and 
Lawrance-Lewis method in Figures (IV-e) and (IV-d) , 
respectively. The Caver method does not have known correla- 
tion but the value of p to give p in the exponential case 
is used. 

As we mentioned earlier, the FXO method has some discon- 
tinuity at the truncation point. But the UXO and TXO methods 
generate relatively smooth continuous distributions. In this 
respect, we say that the FXO is defective and the UXO and 
the TXO are relatively satisfactory. 

The computed correlation from subroutine BIVHST (Bivariate 
histogram) is a little different from given correlation. But 
this can be assiamed as a sampling error. To check this error, 
we simulated 10 times with given correlation p = 0.1. From 
this simulation we get the following results: 



75 



p = mean of computed correlation = 0.105 



VAR[p] = variance of computed correlation = 



o[p] = standard deviation of mean = 0.0064 



0.0004 
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Figure IV-al . Scatter plot for sample of size 2000 from 

mixture-truncation bivariate uniform with 
p = 0.3. Here x^ is fixed at the midpoint 
between the lower and upper bounds. 
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Figure IV-a2. 



Bivariate historgram for sample of 
from mixture- truncation bivariate 
with p = 0.3. Here x is fixed at 
midpoint between the ?ower and the 
bounds . 



size 2000 
uniform 
the 
upper 
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Figure IV-a3. Scatter plot for sample of size 2000 from 

mixture- truncation bivariate uniform with 
p = 0.3. Here x is fixed at the midpoint 
between the lower and the upper bounds. 
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Figure IV-a4 . 



Bivariate histogram for sample of size 2000 
from mixture-truncation bivariate uniform 
with p = -0.3. Here x is fixed at the 
midpoint between the l8wer and the upper 
bounds . 
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Figure IV-bl. Scatter plot for sample of size 2000 from 

mixture-truncation bivariate uniform with 
p = 0.3. Here x is taken to be uniformly 
distributed between the lower and the 
upper bounds . 
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Figure IV-b2. Bivariate histogram for sample of size 2000 

from mixture-truncation bivariate uniform 
with p = 0.3. Here x is taken to be 
uniformly distributed between the lower and 
the upper bounds . 
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Figure IV-b3. Scatter plot for sample of size 2000 from 

mixture-truncation bivariate uniform with 
p = -0.3. Here is taken to be uniformly- 
distributed between the lower and the upper 
bounds 
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Figure IV-b4. 



Bivariate histogram for sample of size 2000 
from mixture-truncation bivariate uniform 
with p = -0.3. Here x is taken to be 
uniformly distributed Between the lower and 
the upper bounds. 
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Figure IV-cl. Scatter plot for sample of size 2000 from 

mixture-truncation bivariate uniform with 
p = 0.3. Here has triangular distribution 
between the lower and the upper bounds . 
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Figure IV-c2. 



Bivariate histogram for sample of size 2000 
from mixture-truncation bivariate uniform 
with p = 0.3. Here x has triangular 
distribution between the lower and the 
upper permissible bounds. 
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Figure IV-c3. Scatter plot for sample of size 2000 from 
mixture- truncation bivariate uniform with 
p = -0.3. Here x has triangular distribu- 
tion between the 2ower and the upper bounds . 
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Figure IV-c4 . 



Bivariate histogram for sample of size 
2000 from mixture-truncation bivariate 
uniform with p = -0.3. Here x has 
triangular distribution betweeS the lower 
and the upper bounds . 
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Figure IV-dl. 



Scatter plot for sample of size 2000 from 
Lawrance-Lewis ' s method with p = 0.3 and 



a 



(2+g) 
3+ (2+6) p 



(| + 1) 



and 6 is uniform over 



2 
3-p 



11 . 
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Figure IV-d2. 



Bivariate histogram for sample of 
from Lawrance— Lewis ' s method with 
and 
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Figure IV-el. 



Scatter plot for sample of size 2000 
from Gaver ' s transformation with p = 0.6. 
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Figure IV-e2. Bivariate histogram for sample of size 

2000 from Gaver ' s transformation with 

p = 0.6. 
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V. BIVARIATE EXPONENTIAL GENERATOR 



Another probability distribution of major interest in 
simulations is the exponential. The cumulative distribution 
function and probability density function for the exponential 
are, respectively. 



The problem of generating exponential deviates reduces to 
one of generating "unit" exponentials, i.e., those with 
X = 1.0, and then multiplying the result by whichever 1/A 
is necessary to give the desired distribution. That is, if 
the random variable E has the exponential (A = 1) distribution 
then. 



F(x) 



1 - e 



- Ax 



X > 0 



and 



f (X) 



A e 



“Ax 



X > 0 



The expected value of the exponential distribution is 



E[X] 



1/A 



and the variance is 



VAR[X] 
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X 




nas the exponential (X = n) distribution. In this section, 
we considered only unit exponentials as marginal distributions 
for bivariate pairs. 

DETERMINATION OF PARAMETERS IN THE EXPONENTIAL MIXTURE- 

TRUNCATION METHOD 

Because X^^ is an exponential (X = 1) truncated to the 
left of Xq and X 2 is also exponential (X = 1) truncated to 
the right of x^, we have 



E[X^] 



0 




1 




VAR[Xj_] 



2 



0 





ElX^] 




00 



F(x) -F(x^) 



X 




b 



1 + X 



o 
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VAR [X^ ] 



X 



» - F(x) - F (x^) 

/ 1-F(x°) 



- y . 



= 1 ; 



jind from definition/ 



1 - a. 



1 1 - ^ “ “2 



= P(x^) 



-X 



= 1 - e 



(V-A-la) 



1 - a . 



2 1 - + 1 “ 0-2 



= 1 - F(x^) 



-X 



= e 



(V-A-lb) 



If we use these formulas in Theorem 2 in Section III, we get 



M = 



"2 2 
— X 
7T O 



(V-A-2) 



3 = — («! - ir,) 

7T 2 1 1 



(V-A-3) 



and 
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p 



(V-A-4) 



77 , O 



For given correlation coefficient p, we can compute as 
a function of from the formula (V-A-4) as 



P 77 



? "1 



-X 



= (1 + -^) (1 - e °) 

X 



and we know that = 1 (1 - a, ) from the formulas (V-A-la) 

2 V 2 J- 

and (V-A-lb) . From this, 



“2 = 1 - 



-X X -X ^ 

= e ° + e °(1 - e 

X 



These and d .2 probabilities, so they have to be greater 

than or equal to zero and less than or equal to one; 



0 < 





< 1 



0 < e ° + e^°(l 



-X 



Oj 2 




< 



1 



(V-A-5a) 

(V-A-5b) 



From these two inequality equations, (V-A-5a) and (V-A-5b) , 
we can find the x^ ranges for given correlation coefficient 
p. To solve these equations, we can divide into two cases. 
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one for positive correlation case and the other for negative 
correlation case. If correlation coefficient p is positive, 
then both equations are always positive. Thus we only need 
to find the ranges which makes (V-A-5a) and (V-A-5b) less 
than 1. From the equation (V-A-5a) , for the case, 

-X 

Oi^ = (1 - e °) (1 + -^) _< 1 

X 

o 

becomes 

p (e - 1) _< x^ 

and let 

_ X 

Z / O T X 

, y 2 = P (e - 1) . 



Because of the first derivatives of y^^ and are always 
positive, we know that these two functions are monotone 
increasing functions. Thus we can find x^ ranges which 
satisfy y^^ — ^2 Newton Raphson method. When using 

the Newton Raphson method, let y = y^^ ~ y2 = 0 find an 

approximate solution, by approximating exponential series, 
which we can use as a starting point. That is. 



X 



X 



y = y^-yj 






3! 



- 1 ) 






X + 

o 



0 
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Then 



X 

o 



(1-^) ± (1- p - ^ 






1/3 p 



Starting with this approximate value in the Newton Raphson 
method, we can find range, say [x^ ,x^ ] , which satisfies 
0 <_ 1. And for the a 2 case, 

-X -X -X - 

= e ° + e °(l-e ^ 1 

^o 

becomes 



X - 

p (e - 1) < X 

— o 

This result is exactly the same as the case, that is, at 

the same range and a 2 satisfy constraints £ 1 and 

Qt~ < 1. Thus we can use x. and x as the lower bound of 
2 - *1 

x^,x , and the upper bound of x ,x . If correlation coeffi- 

cient p is negative, then the equations (V-A-5a) and (V-A-5b) 

are always less than 1. Therefore we need to consider only 

one constraint which makes 0, ot 2 ^ 0. From the 

equation (V-A-5a) , solve the inequality equation 

-X 

0 1 = (1 - e °) (1 + - 2 ^) 
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F 



-X 

since 1 - e ^ is always positive, we see that to satisfy 
the inequality 1 t should be positive, i.e.. 



P 



2 



X 



o 



> 



-1 



equivalently, we have 



X 



o 



> 




In the case, from equation (V-A-5b) , 

-X X “X ^ 

0 = e ° + e °(l-e 

or, equivalently, we have 

_ X ^ 

2 / o 1 ^ 2 

1 -p(e - 1) 

As in the positive correlation case, we can find a starting 
point by approximation to solve this equation by Newton 
Raphson. The result comes out as 



Xg = (/^ + p)/(- j) 

with this starting point we find another bound of x^ which 
satisfies 0 < a-. This becomes the upper bound of x , x , 
and from the case, we have a constraint x^ >_ which 
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becomes the lower bound of x , i.e., x = The lowest 

and highest correlations available for bivariate exponential 

pairs in mixture- truncation method are approximately -0.480 and 

0.647 respectively. By comparison note that the most negative 

correlation available for bivariate exponential pairs with 

identical fixed marginals is 1-^ [Moran (1967)]. Gaver's 

(1972) negatively correlated pair has correlations in the 

range [-0.5,0]. The table (V-1) shows the lower and upper 

bound of X in the mixture-truncation method with identical 
o 

marginal exponential and given correlation. 



Table V-1: The lower bound and upper bound of x 

in the mixture- truncation method with 
identical exponential marginal distributions 



x^ range for each correlation 



p 




X 

u 


P 




X 

u 


0.1 


0.106 


5.832 


1 — 1 

• 

o 

1 


0.317 


1.984 


CM 

• 

O 


0.225 


4.723 


CM 

• 

O 

1 


0.448 


1.439 


0.3 


0.362 


3.990 


1 

o 

• 

U) 


0.548 


1.103 


0.4 


0.527 


3.395 


• 

o 

1 


0.633 


0.855 


0.5 


0.741 


2.842 


-0.45 


0.671 


0.751 


0.6 


1. 082 


2.223 









B. GENERATING PROCEDURE 

The generating procedure of bivariate exponential random 
vector is almost the same as .in the uniform case. As in 
the loniform case, we develop three procedures which are the 
FXO method, the UXO method and the TXO method. As mentioned 
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in Section III, all of these methods are exactly the same 

except in how x^ is chosen from the x^ range The 

FXO method chooses x^ as a fixed point between x^ and x^. 

This procedure makes some discontinuity at the truncation 

point. In this respect, this method is defective. 

The UXO and TXO methods are the same except we use 

different distributions for x . In the UXO method we use 

o 

the uniform distribution and the triangular distribution for 
the TXO method. These methods have more smooth distribution 
than the FXO method. 



Exponential Mixture-Truncation Method 
1. (Initialization) 

i) For given -0.48 < p < 0.64, find x^ and x^ 



2. Define truncation point x^ 

* FXO method 

i) = i(x„ + X ) 

o 2 £ u 

* UXO method 

i) Generate a uniform [0,1] random variable U^^ 



* TXO method 

i) Generate two uniform [0,1] random variables 

V . V 
1 ' 2 

ii) = '^£ ^1 ^2 



where 
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MMPPiilM 














X, 



= (x^ - 



m 






V, 



= (x„ - 



u 






V, 



3. Compute parameter values 

-X 

ir^ = F(x^) = 1 - e ° 

’2 ' ^ * ”l 

“l " 

a. = 1 -a - a,) 

2 7T2 1 

4. Choose type for Y 

i) Generate a uniform [0/1] random variable U 

ii) If U £ 7r^, go to 9 

5. Y is an 

i) Generate an exponential random variable 
ii) If > x^, set Y E^ and go to 6 

iii) Otherwise/ return to 5.i) 

6. Choose type for Z 

i) Set U ^ ((U- 7r^)/(l- 7T^)) 
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ii) If U _< 1 - a^, go to 8 

7. Z is an X 2 

i) Generate an exponential random variable E 2 
ii) If E 2 > set Z E 2 and go to 11 

iii) Otherwise return to 7.i) 

8. Z is an 



i) 


Generate an exponential random variable E 2 


ii) 


If E 2 £ set Z E 2 and go to 11 


iii) 


Otherwise return to 8.i) 



9 . Y is an X;^ 



i) 


Generate an exponential random variable E^ 


ii) 


If E, < x^, set Y ■*- E, and go to 10 
1—0 1 ^ 


iii) 


Otherwise return to 9.i) 



10. Choose type for Z 



i) 


Set U ^ U/tt^ 


ii) 


IfU^a^, go to 8 


iii) 


Otherwise go to 7 



11. Deliver (Y,Z) and go to 4 for the FXO method, or go 
to 2 for the UXO and TXO methods until a sufficient 
number of random vectors are obtained. 

For the exponential case it is possible to give a more effi 
cient algorithm in which and X 2 are generated exactly. 
The algorithm is as follows. 
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Efficient Exponential Mixture-Truncation Method 



1. (Initialization) 

i) For given -0.48 < p < 0.64, find and 

2. Define the truncation point x 

o 

* FXO method 

X = T-(X. + X ) 

o 2 £ u 

* UXO method 

i) Generate a uniform [0,1] random variable 



* TXO method 

i) Generate two uniform [0,1] random variables 



ii) X 



o 





ii) X 



o 




where 



X 



m 




X 



1 




X 



2 




3. Compute parameter values 




1 - e 



o 
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1 



- 7T . 



“l ~ (1 + 

“2 = 1 - “ 1 » 

4. Choose type for Y 

i) Generate a uniform [0,1] random variable U 
ii) If U _< 7 T^, go to 9 

5. Y is an X 2 

i) Generate an exponential random variable 

ii) Set Y X + E, 
o 1 

6 . Choose type for Z 

i) Set U ^ ( (U - 7 T ^)/(1 - 77 ^) ) 

ii) If U _< 1 - tt 2 / go to 8 

7. Z is an X 2 

i) Generate an exponential random variable E 2 

ii) Set Z -e X + E- and go to 11 
o 2 ^ 

8 . Z is an X^ 

i) Generate a uniform [0,1] random variable W 2 
ii) Set Z - ln(1.0 - W 2 *t 7 ^) and go to 11 

9. Y is an X^ 

i) Generate a uniform [0,1] random variable 
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ii) Set Z - ln(1.0 - 

10. Choose type for Z 

i) Set U U/ttj_ 

ii) If U £ go to 8 

iii) Otherwise, go to 7 

11. Deliver (Y,Z) and go to 4 for the FXO method, or go 
to 2 for the UXO and TXO methods until a sufficient 
number of random vectors are obtained. 

Note that to compute and in step 1 of both algorithms 
we use subroutine BOUND which is listed in the appendix and 
the program with effective algorithm also. The scatter plot 
and bivariate histogram of resulting random vectors are shown 
in Section V-D. 

C. REGRESSION OF Z ON Y FOR GIVEN p 

Schmeiser (1979) has used the regression of Z on Y = y 

to fix the parameters in his bivariate gamma distribution. 

Consequently we investigate this for the mixture-truncation 

method case. The regression is different depending on 

whether Y < x^ or Y > x^. We consider two cases here, one 
— o o 

for fixed x^ and the other for x having uniform distribution, 
o o 

For fixed x^, we have 



E[Z[Y = y,y_< x^] 



aj^E[Xj^] + (l-aj^)E[x2l 



1 + X 



o 
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Substituting the value for 



a 



1 




then we have 

E[ZlY = y,y < x^] = 1 + ( 1 + -^) 




And if y > x^, then 

E[Z|Y = y,y > x^] = ( 1 - a 2 ) E [x^ ] + a^E [x^ ] 




Substituting the value for a 2 / 



a 



2 





then we have 



E[Z|Y = y,y > x^] 



1 + — ^ 
2 O 
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Thus the regression is constant over (0,x^) and changes for 
y ^ x^. This is not surprising in light of the joint dis- 
tribution given in Section IV-D. For laniformly distributed 
x^/ the computation is different for different ranges of Y. 
If y _< x^, then we have 



X 

u 

E[z|Y = y] = J E[z|Y = y, X = x^ , 



If 



X 



u 



X, 



X 



I (1 e ° dx^ 



-^u e"'‘° , 

- e + e - p j dx 



X 



2 o 



£ 



-X, -X T -X 

A u . / 1 u 

e - e + p (— e 

X 

u 



1 

— e 



X 

+ In — - X + X . ) 

x„ u ii 



x„ < y < X 



then we have 
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E[Z|Y = y] = / E[Z|Y = y, X=x^, y^x^]f(x^) 



dx 



X 



u 



+ / E[Z|Y = y, X=x^, y_<x^]f(x^) dx^ 

y 



-X 



2e ^ + 2py - p (x^ + x^) + p e ^ ^ 



u 



X 

+ p In — 

y 



If y > X , then we have 
^ u 



X 

u 

E[z|Y=y] = I E[z|Y = y, X = x^ , y>x^]f(x^) dx^ 

X « 



X. 



U TT, -X 

I {1 + — -^) e ° dx 

^ TT^ X O 

2 o 



**X ** X 

e ^ - e '^ + p(y-x) 



By making x^ uniformly distributed over the available range 
of x^ for given correlation, we can get smoother behavior 
for the regression function. 

D. SIMULATION RESULTS 

We will show here the scatter plots and bivariate 
historgram of the mixture-truncation bivariate exponential 
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random vectors with correlation p =0.3 and p = -0.3 by 
the FXO, UXO and TXO methods in Figures (V-a) , (V-b) and (IV-c) , 

respectively. Also we show the results from Gaver's bivariate 
exponential and Marshall and Olkin's in Figures (V-d) and 
(V«e) , respectively. 

In the mixture-truncation bivariate exponential distribu- 
tion, the generated random vectors by the FXO method also 
has discontinuity at truncation point but the other cases 
have relatively smooth distributions. 

The computed correlation printed in the left side of the 
bivariate histogram is a little different from the given 
correlation. But this difference can be assumed as a 
sampling error. 

To check this error we simulated 10 times with given 
correlation p= 0.1. From these simulations we get: 

p = mean of computed correlation = 0.092 

VAR[p] = variance of computed correlation = 0.0008 

aCp) = standard deviation of mean = 0.009. 
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X X 




Figure V-al. Scatter plot for sample of size 2000 from 
mixture- truncation bivariate exponential 
with p = 0.3. Here x^ is fixed at midpoint 
between the lower and°the upper bounds . 
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Figure V-a2, 



Bivariate histogram for sample of size 
2000 from mixture-truncation bivariate 
exponential with p = 0.3. Here x is 
fixed at the midpoint between the^lower 
and the upper bounds . 
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Figure V-a3. Scatter plot for sample of size 2000 from 
mixture-truncation bivariate exponential 
with p = -0.3. Here x is fixed at the 
midpoint between the lower and the upper 
bounds. 
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Figure V-a4. Bivariate histogram for sample of size 
2000 from mixture-truncation bivariate 
exponential with p = -0.3. Here x is 
fixed at the midpoint between the ?ower 
and the upper bounds . 
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Figure V-bl. Scatter plot for sample of size 2000 from 
mixture-truncation bivariate exponential 
with p = 0.3. Here x is taken to be 
uniformly distributed oe tween the lower 
and the upper bounds . 
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Figure V-b2 . 



Bivariate histogram for sample of size 
2000 from mixture-truncation bivariate 
exponential with p = 0.3. Here x is taken 
to be uniformly distributed betwe§n the 
lower and the upper bounds. 
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Figure V-b3. Scatter plot for sample of size 2000 from 
mixture-truncation bivariate exponential 
with p = -0.3. Here x is taken to be 
uniformly distributed Between the lower 
and the upper bounds . 
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Figure V-b4. 



Bivariate histogram for sample of size 2000 
from mixture-truncation bivariate exponential 
with p = -0.3. Here x is taken to be 
uniformly distributed between the lower and 
the upper bounds . 
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Figure V-cl. Scatter plot for sample of size 2000 from 
mixture-truncation bivariate exponential 
with p = 0.3. Here x has triangular 
distribution between ?he lower and the 
upper bounds. 
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Figure V-c2. Bivariate histogram for sample of size 2000 
from mixture-truncation bivariate 
exponential with p = 0.3. Here x has 
triangular distribution between tRe lower 
and the upper bounds . 
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Figure V-c3. Scatter plot for sample of size 2000 from 
mixture-truncation bivariate exponential 
with p = -0.3. Here x has triangular 
distribution between tRe lower and the 
upper bounds. 
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Figure V-c4. Bivariate historgram for sample of size 
2000 from mixture-truncation bi- 
variate with p = -0.3. Here x has 
triangular distribution between the lower 
and the upper bounds. 
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Figure V-dl. Scatter plot for sample of size 2000 from 
Gaver ' s bivariate exponential with 
p = -0.3. 
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Figure V-d2 . Bivariate histogram for sample of size 
2000 from Gaver's bivariate exponential 
with p = -0.3. 
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Figure V-el. Scatter plot for sample of size 2000 from 
Marshall and Olkin's bivariate exponential 
with p = 0.3. 
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Figure V-e2. Bivariate histogram for sample of size 2000 
from Marshall and Olkin's bivariate 
exponential with p = 0.3. 
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VI. BIVARIATE GAMMA GENERATOR 



The gamma distribution with shape parameter r and 
scale parameter X has the density function 



f (X) = 



r (r) 



r-1 -Xx 
X e 



if X > 0 



if X < 0 



(VI-0-1) 



and cumulative distribution function 



F{y) = J 



y X^ 



r (r) 



r-1 -Xx , 

X e dx . 



(VI-0-2a) 



In particular, if the shape parameter r is a positive integer 
then 



F(y) 



1 



r-1 



I 

j=0 



e ^^(Xy) 
jl 




where r (r) is Euler's gamma function 



(VI-0-2b) 



r(r) = / x^~^ e~^ dx . 

0 

If the random variable X has Gamma (r,X) distribution then 
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E[x) 



r 

A 



VAR[x] = ^ 

A^ 

when r = 1, X has the exponential distribution while X, 
suitably scaled, has an asymptotically normal distribution 
as r We note that if X has a Gamma (r,l) distribution 

then has a Gamma (r. A) distribution. So we may set = 1 
as far as the generating algorithm is concerned. The output 
from the generator may then be appropriately scaled. Two 
algorithms for generating bivariate gamma random vectors are 
given in Section II. Since they have complicated computa- 
tions to define the parameter values and use inverse trans- 
formation functions, the mixture-truncation method is 
probably preferable to those methods. To generate bivariate 
random vectors with identical Gamma (r,l) marginal distribu- 
tions and having given correlation p by mixture-truncation 
method, we only need univariate garame (r,l) random variate 
generator. Under the assumption that the univariate gamma 
generator is available we developed this procedure. In fact 
we used Lewis and Robinson's gamma generator as a univariate 
gamma generator. 

A. DETERMINATION OF PARAMETERS IN THE GAMMA 
MIXTURE-TRUNCATION METHOD 

If the marginal distribution is gamma with positive 
integer shape parameter r and scale parameter A = 1, then 
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we know that F(x) becomes cumulative distribution function 

of gamma (r,l). Then the is gamma (r,l) truncated to the 

left of and is also gamma (r,l) truncated to the right 

of X . Thus we have 
o 

X 

EtXil = ° K 



1 

F(x^) r (r) 




X e dx 



1 

(r) 



{r! 




i=0 



rl 

(r-i) 1 



r-i 



]} 



VAR[X^] = 




F(x 

F(Xo 



"l 



2 



X 






r+1 -X , 2 

X e dx - 



1 

(r) 



{ (r+1) 



e 



-X 

o 



r+1 



t I 

i=0 



(r+1) ! 
(r+l-i) I 



r+l-i 

o 



]} 



2 

■ ^1 
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oo 



E[X^] = 



VAR[X2] = 



Let 



A = 



B = 



F(x) - F(x^) 

'2 " ^ ^ 1 - F(x„) 

==0 ° 



— FT — T / e ^ dx 

it2r(r) > 

O 



[e I tTTTTT 1 



^2^ (r) j_£o ° 



2 F(x) - F(x^) 
FU^ 



j ^ d . 



1 f r+1 -X j 

v^J ^ ^ 

o 



1 -==0 _<r±lil ^ r+l-i 2 

H2r(r) (r+l-i)! ° ^ ‘ “2 




r 

(r-i) ! 




i 



X r+1 / , 1 X , .I ■ 

o V (r+1) ! r+l-i 

^£q (r+l-i) ! o 
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and substitute, then we have 



E[X^] = 



r I - A 

TTj^r(r) ' 



VAR[X^] 



= a 



2 ^ (r+1) I - B _ (rl - A) 



TTj^r (r) 



2 2 
TT^^r(r)^ 



ELX^] = 



TT2r(r) ' 



VAR[X2] 



= a. 



B 



ir,r(r) 2„- >2 • 

2 r(r) 



If we use these formulae in Theorem 2 in Section III then 
we have 



and 



(TT2rl - A)^ 





6 = 


ai - TT^ 

/ 

1^2 




Thus 






(VI-A-1) 




P ” 


r 712^ (r (r) ) ^ 


For given p , 


X and 
o 


r. 






^1 = 


2 

p r! (r-1) ! 71^ 772 

2 ’^l 

( 772r ! - A) 


(VI-A-2a) 
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and 



7 T 2 - ^^(1 - a ^) 




(VI-A-2b) 



The value of (VI-A-2a) and (VI-A-2b) are determined with 
given p, and r. Thus for given marginal distribution and 
correlation coefficient p we can find ranges as before. 
For simplicity, we will consider them with marginal distri- 
bution gamma (2,1). From the formula of (IV-0-2b) , 







- X 



o 



e 



-X 

o 



TT 



2 




-X 



+ X 



e 



-X 

o 



and 



-X 



A = 



o / 2 ^ 



2x + 

o 



2 ) 



Then (VI-A-2a) and (VI-A-2b) become 



a 



1 



a 



2 



2 p 



o 




+ TT . 



+ 



2 p 





(VI -A- 3a) 



(VI-A-3b) 



These and are probabilities, so they have to be greater 
than or equal to zero and less than or equal to one; 
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2 



0 



< 



0 



- 2 



2 p 1T^ TT^ 




+ TT. 



< 1 



"2 



^ P ^1 ^2 

. -2x 

4 o 

X e 
o 



< 1 



From these two inequality equations we can find the x^ 

ranges for given correlation coefficient p. If p ^ 0, then 

(x^ and ate always positive. Thus we only need to find the 

X ranges which makes (VI-A-3a) and (VI-A-3b) less than one. 
o 

From equation (VI-A-3a) , the case, and from equation 

(VI -A- 3b) , the 02 case, we have exactly the same inequality 

equation; 



X 

2pe °(1 + X ) 
o 



+ 2p (1 + 



"o' 



(VI-A-4) 



If p £ 0, then and are always less than one. Thus we 
only need to find the x^ ranges which makes (VI-A-3a) and 
(VI-A-3b) greater than zero. From the equation, we have 



v''2T*(1 + Xq) ^ x^^ 



(VI-A-5a) 



where 



P * = I P I 



and from the 02 equation, we have 
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2 



(VI-A-5b) 



X 

/2p * e ^ 



< X 

— o 



+ /^(l + X ) 

o 



where 

p* = 1 P 1 • 

Note that if we let for the left side equation of inequality 
sign and Y 2 for the right side equation of inequality equa- 
tions (VI-A-4) , (VI-A-5a) and (VI-A-5b) , then we know that all 

y^ and Y 2 monotone increasing functions. Thus as in the 

exponential case we can find the range of x^ by the Newton 
Raphson method. For the positive cirrelation case, equation 
(VI-A-4) will give the x^ and x^ and for the negative corre- 
lation case, equation (VI-A-5a) gives x^ and equation 
(VI-A-5b) gives x^. By using subroutine GBOUND which is 
listed in the appendix, we can find x^ and x^. The tables 
(Vl-a) and (Vl-b) show the results from subroutine GBOUND. 
Unfortunately we have limitations on the possible range of 
correlations; we can generate bivariate random vectors with 
correlations in approximately the range -.5 £ p ^ .6 for the 
Gamma (2,1) marginal distribution. In the same way as the 
above, the values x^ and x^ can be computed with increasing 
complexity for integer values of r greater than 2. In 
principle it is also possible, using numerical integration, 
to compute allowable values for non-integer values of r. The 
computation is not as complicated as those required for 
Schmeiser's bivariate gamma. Moreover the mixture-truncation 
method does not require that the user be able to compute 
the inverse gamma distribution function. 
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Table Vl-a: x ranges with Gamma (2,1) 

marginal distribution 
and given p 



p 




X 

u 


P 




X 

u 


0.1 


0.41 


7.38 


-0.1 


0.93 


3.41 


0.2 


0.65 


6.23 


-0.2 


1.18 


2.75 


0.3 


0.89 


5.42 


-0.3 


1.35 


2.34 


0.4 


1.15 


4.74 


-0.4 


1.50 


2.03 


0.5 


1.47 


4.09 


-0.5 


1.62 


1.79 


0.6 


1.95 


3.32 









Table Vl-b: x ranges with Gamma (3/1) 

marginal distribution 
and given p 



P 




X 

u 


P 




X 

u 


0.1 


0.82 


8.99 


-0.1 


1.64 


4.73 


0.2 


1.17 


7.72 


-0.2 


1.37 


3.98 


0.3 


1.49 


6.82 


-0.3 


2.21 


3.5 


0.4 


1.83 


6.05 


-0.4 


2.4 


3.15 


0.5 


2.24 


5.31 


-0.5 


2.56 


2.89 


0.6 


2.84 


4.42 
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B. GENERATING PROCEDURE 



The generating procedure of bivariate gamma random vec- 
tors is almost the same as the uniform and exponential cases. 
We also developed here three methods which are the FXO 
method, the UXO method and the TXO method. As in the uniform 

and exponential cases, these three methods are the same 

except in the choice of x from the x range [x.,x ]. The 
FXO method chooses x^ as a fixed point from the x^ range and 
uses it during the entire routine while the UXO and TXO 
methods choose another x^ in every routine. From experience, 
the midpoint of the x^ range gives the best result of the 

FXO method. In the UXO method and the TXO method, we assume 

that x^ has the uniform distribution and triangular distri- 
bution, respectively, over [x^,x^]. This assumption removes 
some discontinuity which occurs at the truncation point in 
the FXO method. 

Gamma Mixture-Truncation Method 
(for integer shape parameters) 

1. (Initialization) 

i) For given allowable correlation coefficient p. 



find X. and x 

Z ■ 



u 



2. Define truncation point x 



o 



* FXO method 



i) X 



o 
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* UXO method 



i) Generate a uniform [0,1] random variable 
ii) 

* TXO method 

i) Generate two uniform [0,1] random variables, 

Vi, V2 



ii) X = X + X- + x„ 

o z 1 2 



where 



X = (x + X ) /2 . 0 , 

m z u ' 



=^1 = ' 



X- = (x - X ) *V- 
2 u m 2 * 



3. Compute parameter values 



-X 

r-1 e ° X ^ 

F(x ) =1-1 ° 



j=0 



I f 



= 1 - 



p r ! (r-1) ! 7T^ Ti 2 






2 * *1 ' 



’1 

1 - “i’ - 



where 
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A 



e 



r 



r-i 



-X 



r 

I 

i=0 



(r-i) 



X 



4. Choose type for Y 

i) Generate a uniform [0/1] random variable U 

ii) If U ^ TTj^/ go to 9 

5. Y is an X 2 

i) Generate a gamma (r,l) random variable G^ 

ii) If Gt > X / set Y G, and go to 6 
1 o 1 ^ 

iii) Otherwise return to 5.i) 

6. Choose type for Z 

i) Set U ^ ( (U - TT^)/(1 - 77j^) ) 

ii) If U _< 1 - a 2 # go to 8 

7. Z is an X 2 

i) Generate a gamma (r,l) random variable G 2 

ii) If G 2 > x^/ set Z G 2 and go to 11 

iii) Otherwise return to 7.i) 

8. Z is an X^ 

i) Generate a gamma (r,l) random variable G 2 

ii) If G_ ^ X / set Z G„ and go to 11 

iii) Otherwise return to 8.i) 
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9 . Y is an 

i) Generate a gamma (r,l) random variable 

ii) If G^ _< x^, set Y -f- G^^ and go to 10 

iii) Otherwise return to 9.i) 

10 . Choose type for Z 

i) Set U u/tt^ 

ii) If U £ a^, go to 8 

iii) Otherwise go to 7 

11. Deliver (Y,Z) and go to 4 for the FXO method, or go 
to 2 for the UXO and TXO methods until a sufficient 
number of random vectors are obtained. 

For the non-integer shape parameter case, step 1 and step 
3 need more complicated computation procedures but here we 
are only concerned for integer shape parameter cases. In 
step 1 we used the subroutine GBOUND which is listed in the 
appendix. The scatter plots and bivariate histograms of the 
resulting random vectors from this algorithm are shown in 
the next section. 

C. SIMULATION RESULTS 

The scatter plot and bivariate histograms for random 
vectors of the mixture-truncation bivariate gamma variables 
with marginal gamma (2,1) distribution and given correlation 
(0.3 and -0.3) are given here. As in the uniform and 
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exponential cases, the UXO, TXO and FXO methods are used 
and are shown in Figures (Vl-a) , (Vl-b) and (VI-c) , respectively. 

The FXO method has some discontinuity at the truncation 
point, but the UXO and TXO methods appear to give a relatively 
smooth continuous distribution. The computed correlation in 
subroutine BIVHST (Bivariate histogram) is a little differ- 
ent from the given correlation. But we can assume this 
difference is a result of sampling error. 
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X 



X 



X 




Figure Vl-al. Scatter plot for sample of size 2000 

from mixture-truncation bivariate Gamma 

distribution with p = 0.3. Here x has 

uniform distribution over (x„,x ). ° 
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Figure VI-a2. 



Bivariate histogram for sample of size 2000 
from mixture-truncation bivariate Gamma 
distribution with p = 0.3. 
uniform distribution over 



Here x has 
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Figure VI-a3. Scatter plot for sample of size 2000 from 

mixture-truncation bivariate Gamma 
distribution with p = -0.3. Here x has 
uniform distribution over (x ,x ). ° 
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Figure Vl-bl. Scatter plot for sample of size 2000 

from mixture truncation bivariate Gamma 
distribution with p = 0.3. Here x has 
triangular distribution over (Xjj^,x°). 
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Figure VI-b2. 



Bivariate histogram for sample of size 
2000 from mixture-truncation bivariate 
Gamma distribution with p = 0.3. Here 
X has triangular distribution over 
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Figure VI -b3. Scatter plot for sample of size 2000 from 

mixture- truncation bivariate Gamma 
distribution with p = -0.3. Here x has 
triangular distribution over (x 
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Bivariate histogram for sample of size 
2000 from mixture- truncation bivariate 
Gamma distribution with p = -0.3. Here 
X has triangular distribution over 
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Figure VI-cl. Scatter plot for sample of size 2000 from 

mixture-truncation bivariate Gamma 
distribution with p = 0.3. Here x^ is 

fixed at the midpoint of x and x . 
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Figure VI-c2. 



Bivariate histogram for sample of size 

2000 from mixture-truncation bivariate 

Gamma distribution with p = 0.3. Here 

X is fixed at the midpoint of x^ and 
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Figure VI-c3. 



Scatter plot for sample of size 2000 
mixture truncation bivariate Gamma 
distribution with p = -0.3. Here 
fixed at the midpoint of x and 
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Figure VI~c4. 



Bivariate histogram for sample of size 
2000 from mixture- truncation bivariate 
Gamma distribution with p = -0.3. Here 
is fixed at the midpoint of x and x 
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VII. CONCLUSION 



The mixture-truncation method is a general method which 
can generate bivariate random vectors having any theoretical 
marginal distribution and allowable correlation. The 
generating procedure is very simple and doesn't need much 
computation for defining parameter values. In this respect, 
the mixture-truncation method is a very attractive method for 
generating bivariate random vectors. A price is paid for 
this simplicity and generality in that the Frechet bounds of 
correlation for the bivariate distributions specified by the 
marginal distribution given by Moran (1967) are not always 
attained. Also there is some discontinuity in the bivariate 
distribution. However this discontinuity can be decreased 
by giving some distribution to the truncation point over 
its range for given p . Thus the mixture-truncation method 
is very attractive for simulation studies involving only 
partly specified dependency structures. The mixture- 
truncation method may be extended to generate bivariate 
random vectors having negative values. Another extension 
may be made to use grade correlation or rank correlation 
which are invariant under transformation instead of using 
the product moment correlation. 
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APPENDIX: PROGRAM LISTINGS 



THE PROGRAM FOR GENERATING 3IVARIATE RAOQM V=CTCRS(Y,Z) 
HAVING IDENTICAL MARGINAL DISTRIBUTION BY MI X'URE-'^R'JNCA 
C-TICN METHOD. 

THE CORRELATION LIMItc ARE AS FOLLOWS FOR GIVEN MARGINAL 

distribution. 

1. LNIFORM; -0.75 < OR < 0.75 

2. EXPONENTIAL : -0.48 < CP < 0.64 
i.GAMMA(2 ,i ) ; -.55 < OR < 0.64 

GENERATING PROCEDURE 

CALL BUNP(CR,NV,Y,Z,TS) 

CALL B5XP(CRtNV,Y, Z,IS) 

CALL BGAM2(CRtNV, Y, ZtlS) 

WHERE 

OR; GIVEN CORRELATION. 

NV; NUMBER OF RANDOM VECTORS BE GENERATING. 

Y; THE FIRST ELEMENT OF VECTOR WITH DIMENSION NV. 

Z; THE SECOND ELEMENT OF VECTOR WITH DIMENSION MV. 

IS; INITIAL SEED FQR UNIVARIATE GENERATOR. 



SUBROUTINE BUNF ( CR » NV ,Y , Z , lU ) 

C 

C 

CIMENSION Y(NV),Z( NV) 
r 

IF(CR .GT.0.75 .OR. CR.LT.-0.75) GO TG 70 

1= (CR.LT.0.0 ) GO ^0 5 

XL = C.5-( ( (9.0-12.0*CR)*=«‘0.5)/6.0) 

XU=1.0-XL 
GC TO 10 

5 XL = ( (-1 .0*CR) /3 .0 )*=*0.5 
XU=1.0-XL 
10 CO 100 1=1, NV 

Call ranoom(il,u,i) 

XO=XL + ( XU— XL )^U 

CALL UPJOBICR ,X0,AP1, AP2, ?I E 1 , P I E2 , 8 ETA ) 

C CHOOSE TYPE FOR Y 

CALL RANDOM( IU,U ) 

IF(L .LE.PIEl ) GO TO 50 
L=(U-PIEl)/( 1.0-PIEl) 

AP=1.C-AP2 
C Y PROM X2 

20 CALL RANDOM! !L,W,: ) 

Y ( I ) =XO+W*( 1 .C-XO ) 

C CHCJSE type for z 

IF(U.LE.AP) GC TO 40 
C...Z FROM X2 

30 CALL RANDOM! IL,W, 1) 

Z!I )=XO+W*!1.0-X0) 

GO TO 100 
C...Z FROM XI 

40 CALL RANDOM! IL,W, 1) 

Z!I )=XO>»=W 
GO TO 100 
C....Y FROM XI 
50 U=U/PI£1 

60 CALL RANDOM! IU,W ,1 ) 

V! I ) = XO=(‘W 

IF!U.LE.AP1) GO TO 40 
GO TO 30 
100 CONTINUE 
RETLRN 

70 WRITE!6,80)CR 

30 F0RMAT!5X, 'CHECK CORRELATION LIMIT, CO NOT ALLOW GIVEN, 
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CC3PRELATICN IN THIS ORCGRAM FOR CR=»tP7.3) 

STOP 

ENC 

C 

C 

SUBROUTINE UP JOB ( C R , XO , AP 1 , AP 2 » P I E 1 , P I E2 » BE^ A ) 
C UNIFOpy 

C FIND INITIAL VALUE 

fioi=<CR + 3.0^XC**2)/(3.0>»X0) 

AP2 = 1.0-( 1.0-APl)’*' (XQ/ (l.O-XO) ) 

PIE1=( 1.0-AF2)/( 1.0-AP1+1.0-AP2) 
FIE2=(1.0-AP1)/(1.0-AP141 .0-AP2 ) 

8FTA=AP1-( 1 .0-AP2 ) 

RETURN 
ENC 



SUBROUTINE BEXP(CR» NV,Y,Z,IS» 

C 

C 

DIMENSION Y(NV) tZ(NV) 

C 

IX=IS*0.8 
IU=IS*0. 3 + 123456 

IF (CR .GT.0.6A .OR. CR.LT.-0.48) GO TO 70 
call B0UNC(CR tXSl tXS2»XL»XU) 

CO IOC 1=1, NV 
CALL RANDOM! IU,U, 1 ) 

XO = XL + ( XU-XL) ^>0 

call EPJ0E(CR ,XO,A01,AP2,PIE1,PIE2,EETA) 

C CFCOSE TYPE FOR Y 

CALL RANDOM! IU,U, 1 ) 

IF!L.LE.P!E1 ) GO TO 50 
L= !L-°IE1 )/! 1 .0-PIEl) 

AP=1.0-AP2 

C Y FR0MX2 

2C CALL EXP0N!IX,X,1) 

V! I )=X + XO 

C CFCOSE TYPE FOR Z 

IF!L.LE.AP) GC to 40 

C Z from X2 

30 CALL EXP0N!IX,X,1) 

Z!I )=X+XO 
GC TC 100 

C Z FROM XI 

40 CALL RANOOM!!L,U, 1) 

2! I )=-l .OYALOG!! .0-U+PIEl ) 

GO 10 100 

C Y FROM XI 

50 U=U/PI51 

60 CALL RANDOM! IL,U,1 ) 

Y !I )=-1.0YAL0G! 1.0-UYPIEl) 

IF!L.LE.AP1 ) GC TO 40 
GO TC 30 
100 CONTINUE 
RETURN 

70 V<RITE!6,8C)CR 

80 F0RMAT!5X , 'CHECK CORRELATION LIMIT, DC NOT ALLCW GIVEN, 
CCOPPELATICN IN THIS PROGRAM FCR CP=',P7.3) 

STCF 

END 

C 

C 

SUBROUTINE EPJOB!CR,XO, API , AP 2 , P I El , P I E2, BETA ) 

C EXPONENTIAL 

PIEl = i.0-!1.0/EXP!XO) 

PIE2=1.0-PI El 
AP1=PI£1*!1 .0 + !CR/XCYY2 ) ) 

AP2=1.0-! ! PIE 1/PI E2)Y( 1.0-APl ) ) 
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B=Ta=APl-(1.0-AP2) 

RETURN 

ENC 



SUBROUTINE BOUND ( RO ,XS1,XS2,XL,XU) 
IF(RC.LT.O.O) GO TO 10 
AA=(1.0-RC-( ( R0**2) *(5.0/12.0 ) ) )**0.5 
XS1=( (1.0-0.5*R0)-AA)/(R0/3.0) 

XS2 = ( (1 . 0-0.5 ♦ R0)+A^) /(RO/3.0 ) 

CALL RPSON( XS1»XX,RC) 

X L= X X 

CALL RPS0N(X52,XX,RC) 

XU=XX 

RETURN 

10 AR=A0S(RO) 

XL=AR**0. 5 
XS1=C.0 

XS2=(( AR**0.5 )-AR)/ (AR*0.5) 

CALL RPSON( XS2»XX,RG) 

XU = XX 

RETURN 

END 



SUBROUTINE RPSQN ( XS , XX , PO) 

EE=C.0001 
Y1 = XS 
N = 1 

10 IF( RG.LT.0.0) GO TO 20 
CALL FUNP(Y1»FV,DV»RC) 

GO TO 25 

20 CALL FUNN(Y1, FVtOVtRO) 

25 Y2=Y1-(FV/DV) 

IF( (ABS( Y1-Y2 )) .LE . EE) GO TO 30 
Y1 = Y2 
GO TG 10 
30 XX=Y2 
RETURN 
END 
C 
C 

SUBROUTINE FUNP( Y 1 » FV . OV » RG ) 

FV=(Y1**2 )-(RG*(EXP(Yl)-1.0) ) 

CV=(2.0*Y1 )-( RQ*EXP(Y1 ) ) 

RETURN 

END 

C 

C 

SUBROUTINE FUNN( Y1 , F V , DV , RO ) 

RR= A E S ( RO ) 

FV = (Y1**2 )-PR*( (EXP (Y1 )-1.0 )**2 ) 

CV=2.0*Yl-( 2. C*RR*EXP( 2 .0*Y1) ) + ( 2 .0*RR *EXP ( Y 1 ) ) 
RETURN 
END 



SUBROUTINE BGAM2 ( C R t NV ,Y , Z , I S ) 

C 

C 

DIMENSION Y (NV) ,Z(NV) 

C 

IX=IS*0.8 
IU=IS*0. 3 + 123456 

IF(CR.GT.0.64 .OR. CR.LT.-0.55) GO TO 70 
CALL GP0UND(CRtXSl.XS2»XL»XU) 

CG 100 I=ltNV 
CALL RANOCM( IL»U,1) 

XQ=XL+( XU-XL) 
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CALL GPJ08(CP ,X0,AP1,AP2,PI51,PI E2, EETA) 

C CHCaSE TYPE FOR Y 

CALL RANOCM< ILtUtl ) 

IF(L.LE.PIEl) GO TO 50 
U=(L-PIE1)/(1.0-PIE1) 

AP = 1 .0-AP2 
C Y FRCNX2 

20 CALL GAMA(2.C,IX,X,1) 

IF(X.LE.XC) GC TO 2C 
Y( I )=X 

C CFCOSE TYPE FOR Z 

IF(U.LE.AP) GC TO AO 

C Z FROV X2 

30 CALL GAMA(2.0,IX,X,1) 

IF(X.LE.XC) GC TO 3C 
Z( I ) = X 
GO TO 100 

C Z FRCP XI 

AO CALL GAMA(2.0,IX,X,1) 

IF(X.GT.XC) GC TO 4C 
2( I }=X 
GO TO 100 

C Y FROP XI 

50 L=U/PIE1 

60 CALL GAMA (2.0,IX,X,1 ) 

IF(X.GT.XO) GC TO 6C 
Y( I ) = X 

IFa.LE.APl ) GO TO AO 
C-G TC 30 
100 CONTINUE 
RETURN 

70 V.RITE(6t80)CR 

80 FaRVAT(5X » 'CHECK CORRELATION LIMIT, CO NOT ALLCVi GIVEN, 
CCORRELATICN IN THIS PROGRAM FOR CR=',F7.3) 

STOP 

END 

C 

C 

SUBROUTINE GP JOB ( CR , XO , AP 1 , AP 2, P I E 1 , P I E2 , BE ta ) 

C GANMA(2,X) 

EX=1.0/EXP(X0) 

FIE1=1 .0-EX-XC+EX 
FIE2=1.0-PIE1 

AD=(X0*=*A)*(1.0/EXP(2.0*X0) ) 

APl=PIEl + ( < 2.C*CR*PIE1*PIE2=*‘*2)/ AC ) 

AP2=1.0-( (PIE1/PIF2)*( 1 .0-APl ) ) 

BETA=AP1-(1.0-AP2) 

RETURN 

END 

C 

C 

SUBROUTINE GB CUND ( RC ,X S 1 ,XS 2 , XU , XU' ) 

IF(RO.LT.O.O) GO TO 10 

Al=( (A.0*R0)/3,0) +SQRT( ((A.O*RO«*2 ) / 9 .0 ) + ( A . 0 *R0 ) ) 
A2=2.0*( 1 .0-( RO/3.0)) 

XS1=A1/A2 
XS1=XS1+1 .0 
>S2=XS1+1 C.O 
CALL GRPS0N(XS1,XX,R0) 

XL = XX 

CALL GRPSONC XS2, XX,PC) 

>U = XX 
RETURN 

10 AR=2.0*ABS(R0) 

XL=(SGRT( AR) + SQRT( AR+SQ RT ( AR ) «A.O ) )/2.0 
XS1=C.0 

A1=SQRT( (SQRT(AR)/6.0)+(R0/9.0) )- ( SCRT( AR ) /6. C ) 

A2=SGRT(AR)/12.0 

XS2=A1/A2 

CALL GRPS0N(XS2,XX,P0) 

XU = XX 
RETURN 
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C 

C 

SU8FCUTINE GR FSON ( X S , XX ,R0 ) 

EE=C.001 

V1 = XS 

10 IF (RC.LT.0.0) GO TO 20 
CALL GFUNP( Y1 ,FV,0V,R0) 

GC TG 25 

20 CALL GFUNNtYl ,FV, DV,PC) 

25 Y2=V1-(FV/DV) 

IF( (ABS(Y1-Y2)) .LE.EE) GO TO 30 
Y1 = V2 
GO TQ 10 
30 XX=Y2 
RETURN 
END 
C 
C 

SUEPOUTINE GFLNP ( Y1 ,FV,DV,RO) 

Fl= < Yl**4) + (2 .0*RO* ( 1.0+Y1)’**2.0 ) 
F2 = 2.0*R0*EXP(Y1 )*( 1 .0 + Yl ) 
FV=F1-F2 

CF1 = (4.0*Y1**3) + (4.C=»R0*( 1.0+Yl) ) 
CF2=2.0*RG*EXP(Y1)*(2.0+Y1) 
CV=CF1-DF2 
RETURN 
END 
C 
C 

SUEFCUTINE GFUNN ( Y1 , FV» DV , RO ) 

TR = SQRT{2.0*A3S(R0) ) 
F1=(Y1*’«'2.0)+(TR:^( 1.0+Yl) J 
F2 = TR*EXP(Y1 ) 

FV=F1-F2 
CF1=(2.0*Y1 )+TR 
CF2=TR:«t£xP( Y1 ) 

CV=CF1-0F2 
PETLRN 
END 



SUEPOUTINE BIVHST (X» Y, N, WORK) 

C 

IMPLICIT REAL’^a (D) 

INTEGER+2 CELL, MAX, SCALE, ONE, CNE5 
INTEGER RUNS, WN, SN, WEIGHT, U, FLAG 

PEAL D, DELTA, DELTAX 

C 

DIMENSION X(N), Y(N), WCRK(N) 

DIMENSION XM0V(6), YM0M(6), OWCRKIE), CELL(32,32) 
DIMENSION KEY(16) ,XLA8EL(8) 

DATA NOUT/6/, ONE/1/, 0NE5/15/ 

C 

IF(N .GT. 15) GO TO 20 
WRITEINOUT, 10) N 

10 FORMAT! •1***BIVHST4'** TOO FEW SAMPLE POINTS. N =',IS) 
RETURN 
C 

20 AN = N 

AN2 = AN * AN 
NM = N - 1 

C CRDER X-SAMPLE AND FIND RANGES 

CALL SORTON 0,1, N,Y) 

XR ANGE = X( N) - X( 1 ) 

VMA> = Y( 1) 

YM IN = Y( 1) 

CO 30 1=2, N 

IF(Y(I) .GT. YMAX) YMAX = Y(I) 
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IF(Y(I) ,LT. YMIN) YMIN = Y(I) 

30 CONTINUE 

YR^NGE = YMAX - YMIN 

1P(XRANGE * YRANGE .GT. 0.) GO TO 5C 
WRITE(N0UT»40) 

40 FORMAT! '1***BIVHST*** X OR Y SAMPLE HAS ZERO RANGE') 
RETURN 
C 

C 2ERC OUT CATA ARRAY 

5C DO 7C 1=1 »32 
CO 6C J=l»32 
CELL(I,J) = 0 
60 CONTINUE 
70 CONTINUE 

C FINE SCALE FACTORS 

AX = 31.999 / XRANGE 
ex = 1.0001 - AX * X(1 ) 

C5LTAX = XRANGE / 32. 

AY = 31.999 / YRANGE 
EY = 1.0001 - AY * YMIN 
CELTA = YRANGE / 32. 

MAX = 0 

C ACCUMULATE CELL COUNTS 

DO 100 1=1, N 
IX = AX * X( I ) + BX 

!Y = AY * Y(I ) + BY 

lY = 33 - lY 

CELLdX, lY) = CElLdX, lY) + ONE 
IFICELLdX, lY) .GT. MAX) MAX = CELLdX, lY) 

100 CONTINUE 

C SCALE FACTCR FCR COUNTS 

SCALE = (MAX - ONE) / ONES + ONE 

C FINC SAMPLE MOMENTS 

QUORK(l) = O.C 

CALL ACCUM(X,N,DWORK) 

CALL MQMENT(OWORK, XMOM) 

OWCPK(l) = O.C 

CALL ACCUM(Y,N,OWORK) 

CALL MOMENT (DUORK. YMOM) 

C FINC COVARIANCE 

CSUM =0.0 
DO 110 1=1, N 

CSUM = OSUM + (Xd) - XMOMd)) * (Y(I) - YMOM(D) 

110 CONTINUE 

COVAR = DSUM / AN 
XDEV = SORT( XMOM( 2) ) 

YOEV = SQRT(YMQM(2) ) 

RHO = COVAR / (XOEV * YDEV) 

C KENCALL'S TAU COEFFICIENT 

T = C.O 

IF(N .GT. 500) GO TC 119 
CO 118 1=1, NM 
DO 114 J=I,N 
2 = Y(J) - Yd) 

IF(Z .GT. 0.) T = T + 2.0 
IF(Z .LT. 0.) T = T - 2.0 
114 CONTINUE 

118 CONTINUE 

T = T / (aN2 - AN) 

C SET UP WORK AREA AND ORDER Y SAMPLE 

119 CONTINUE 
VI = 1.0 

CO 120 1=1, N 
WORK (I) = WI 
VI = WI + 1.0 

120 CONTINUE 

CALL SORTON(Y, 1, N, WORK) 

C FINC RANK'CRDER S'^ATISTICS 

VN = 0 
SN = C 
RUNS = 0 

IF( Y( 1) .LT. Ml ) ) RUNS = 1 
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130 



140 



AN 



N) GO TO ZOO 



150 



SN = 
IX = 



160 



^ 170 



180 



190 



FLAG = 0 
C = 0.0 
! = 1 
IX = 1 
lY = 1 

Z = ABS(FL0AT(!X-TY)) / 

IF(C .LT. Z) C = Z 
IP(X(IX) - Y(IY)) 150, 160, 140 
•CURRENT VALUE IN POCLEO SAMPLE IS A Y 
RUNS = RUNS + flag 
FLAG = 0 
lY = lY 1 
1=1 + 1 
IF(IY .GT 
GO 1C 130 

•CURRENT VALUE IN PCCLED SAMPLE IS AN X 
RUNS = RUNS + 1 - FLAG 
FLAG = 1 
I^N = WN + I 

SN + WEIGHT (I ,N) 

IX + 1 

IFdX .GT. N) GO TO ZZO 
1 = 1 + 1 
GO TO 130 

•X AND Y VALUES ARE TIED. USE MIDRANK CORREC^ICN 
BASE = X( IX) 

N = 1 

J1 = I + I + 1 

J2 = WEIGHT(I,N) + W5IGHT(I+1 ,N) 

K = I 
1 = 1+1 

■GET TIED X VALUES 
!X = IX + 1 

!F( IX .GT. N) GO TQ 180 
IF(X(IX) .NE. BASE) GO TO 180 
I- = N + 1 
1 = 1 + 1 
J1 = J1 + 

J2 = J2 + 

GO TO 170 

•FIND TIED Y VALUES 
lY = lY + 1 

IF ( lY .GT. N) GO TO 190 
IF(Y(IY) .NE. BASE) GO TO 190 
1 = 1 + 1 
J1 = J1 + I 
J2 = J2 + WEIGHT! I ,N) 

GO TO 180 

•WEIGHTED AVERAGE OF TIED VALUES 
L = I - K + 1 

Jl) / FLCAT(L) + 0.5 



WEIGHT! ! ,N) 



J2) / FLOAT!L) + 0.5 



! L-M) ) / =LOAT!L) 



2 = FLOAT !M ♦ 

WN = WN + Z 
Z = FLQAT!M * 

SN = SN + Z 

AO FCC CORRECTICN FCR RUNS 

RUNS = PUNS + FLOAT! 2 * M * 

1 = 1+1 

IF!IX .GT. N) GO TO 220 
IF!IY .LE. N) GO TO 130 

Y CESERVATICNS EXHAUSTED 

200 RUNS = RUNS + 1 
N2 = N + N 
DO 210 1=1, N2 
WN = WN + I 
SN = SN + WEIGHT! I,N) 

210 CONTINUE 

FINC MANN-WHITNEY STATISTIC 

220 L = WN - !N * !N + 1) ) / 2 

SPEARMAN'S RANK CORRELATION COEFFICIENT 

CSUM = 0. 

WI = 1.0 
DO 230 1=1, N 
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CSUM = DSUM + (WI - V^CRK(D) 2 
V«I = WI + 1.0 
230 CONTINUE 

R = 1. - 6. * DSUM / (AN * (AN2 - 1.)) 

C FINC sample medians 

IF(N0D(N,2) .EQ. 0) GO TO 240 
M = 1 + N/2 
XMEC = X(M) 

TMEC = Y( V) 

GO TO 250 
240 t- = N/2 

XMEC = 0. 5 * (X( M) + X( M-f-1) ) 

YMEC = 0.5 * (Y(M) + Y(M+1) ) 

C FIND NORMALIZED STATISTICS 

250 PNCPM = (RUNS - N - 1) * SQRT((AN + AN - 1.) / (AN2 - 
CAN) ) 

FAC = SQRT(12. / ( AN2 * (AN + AN + 1.))) 

UNORM = (U - C.5 * AN2) * FAC 

VNOPM = (WN - AN2 - 0.5 * AN) * FAC 

SNCPM = ( SN - AN2 - 0.5 * AN) * FAC 

C REOROER Y SAMPLE 

CALL SORTON (V.ORK, 1, N, Y) 

VRITE OUTPUT 

FEACING 

WRITE(NOUT, 300) N 

300 FQR^^AT( '1' .20X, 'BIVARIATE SAMPLE SUMMARY' ,10X SAMPLE 
C SIZE =', 19// ) 

V»RITE(N0UT,31C) 

310 F0RvaT( 13X, '+ — ' ,7('**,7('-')),' * +• ) 

C 

C ECD'r OF PLOT 

CALL OUTPUT(CELL( 1 , 1), SCALE » 1, 0.) 

CALL REPEAT 

CALL OUTPL'T(CELL( 1 , 2), SCALE, 2, VMAX - DELTA) 

UR IT£(N0UT,32C) 

320 FORMAT( '+' ,93X, 'MEASURES OF ASSOCIATION') 

CALL REPEAT 

CALL OUTPUT(CELL( 1 , 3), SCALE, 3, C.) 

URITE(NCUT,330) COVAR 

330 FORMAT! ' + ' , 83X , ' C OV ARI ANCE ' ,2 5X , IP E 1 4 .6 ) 

CALL REPEAT 
VRITE(N0UT,34C) RHC 

340 FORMAT! '+• ,83X, 'CORRELATION COEF F IC I ENT ' , 13XFS . 6 ) 

CALL CUTDUT(CELL(1 , 4), SCALE, 4, 0.) 

URITE(NOUT, 350 R 

350 FOPMATC + ' ,83X,' SPEARMAN"S RANK CCRRELA^ICN CCEF. 
C,F9.6) 

CALL REPEAT 

IF ( N .LE. 500 ) WPITE(N0UT,360) T 
360 FORMAT! 83X, 'KENDALL" S TAU C CE FF I C I ENT ' , 1 IX F9 .6 ) 
CALL 0UTPUT(CELL( 1 , 5), SCALE, 5, C.) 

CALL REPEAT 

CALL OUTPUT(CELL( 1 , 6), SCALE, 6, YMAX - 5 . =♦ DELTA) 
URITE(N0UT,37C) 

370 FORMAT! ' + ' .95X, 'TESTS FOR EQU I D IS T P IBUT ION ' ) 

CALL REPEAT 

CALL OUTPUT (CELL!! , 7), SCALE, 7, 0.) 

URITE!N0UT,38C) 

380 FORMAT! '+', 11 IX, 'TEST' , 7X, ' NORMAL I ZED* ) 

CALL REPEAT 
URnE!NOUT, 390 

390 FORMAT! ' + ', 109X, * STATISTIC* , 5X , * ST AT I ST I C * ) 

CALL OUTPUT!CELL! 1 , 8), SCALE, 8, 0.) 

URITE!NOUT,40C) 0 

400 FORMAT! ' + ' ,83X,' KOLMCGCPOV-SMIRNOV TEST *,F10.6) 

CALL REPEAT 

WR ITE!N0UT,41C) RUNS, RNORM 

410 FORMAT! '+',83X,'WALO-WCLFOWITZ RUNS TEST ' , 1 10 , 5XF9 . 5 ) 
CALL OUTPUT!CELL( 1 , 9), SCALE, 9, 0.) 

URITE!N0UT,42C) U, UNORM 
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420 FORh'AK • + • »83>,’MAN^-WHITNEY TEST • , 7X 1 10 , 5X F9 .5 ) 



430 



440 



450 



460 



470 



480 



490 



500 



510 



520 



530 



SCALEtllt 0.) 



CALL repeat 
V»RITE(N0UT ,430 WN, WNCRM 

FQRyAT( •+*, 8 3X,' WILCOX SOM TE S T» , 1 OX 1 10 , 5X >=9 . 5 ) 

CALL QUTPUT(CELL(1 ,10) , SCALE, 10, YNAX - 9. =* DELTA) 
k»RITE(N0UT,44C) SN, SMCRM 

F□R^<AT( • + • ,83>, 'STEGEL-TUKEY TEST* , 7X1 10, 5XF9 .5 ) 

CALL REPEAT 

CALL OUTPL'T(CELL( 1 ,11 ) 

CALL REPEAT 
V<RITE(N0UT,45C) 

FORMAT! • + • ,92 UNIVARIATE STATISTICS* ) 

CALL OUTPUT(CELL( 1,12), SCALE, 12, C.) 

CALL REPEAT 
V»RITE(N0UT,46C) 

FCPMAT( •+* ,95X, *X S AMPL E* , 9X, * Y SAMPLE') 

CALL OUTPUT(CELL( 1 ,13) , SCALE, 13, 0.) 

V»RITE(N0UT,47C) XMOM(l), YMQM(l) 

FORMAT! *-•- * ,83X,*MEAN*, 5X1P£14.6,3XE14.6 ) 

CALL REPEAT 

WRITE!N0UT,48C) XMEC, YMEO 

FORMAT! * + • ,83X,* MIEOIAN* ,3X1PE 14. 6 , 3X E 14. 6 ) 

CALL OUTPUT!CELL! 1 , 14) , SCALE, 14, YMAX -13. =» DELTA) 
CALL REPEAT 

kRITE!N0UT,49C) XM0M!2), YM0M!2) 

FORMAT! *+* ,83X,*VARIANCE* ,1X1PE14.6,3XE14.6) 
call OUTPUTICELL! 1 ,15) , SCALE, 15, C.) 
t»RITE(NOUT,50C) XOEV, YOEV 

FORMAT! *+ * ,83X, *STO DE V * , 2X IP El 4. 6 , 3 XE14. 6 ) 

call repeat 

^«RITE(NOUT, 510 XRANGE, YRANGE 
FCRMAT!*+*,83X,*RANGE*,4X1PE14.6,2XE14.6) 

CALL 0UTPUT!CELL( 1 ,16) , SCALE, 16, 0.) 

CALL REPEAT 

V»RITE!N0UT, 52 0 XMCM!4), YM0M!4) 

FORMAT! »+• ,83X,* SKEWNESS* ,1 X1PE14 .6 ,3XE14.6 ) 

CALL OUTPUT!CELL( 1 ,17) , SCALE, 17, 0.) 



CALL OUTPUT!CELL( 1 ,17) , SCALE, 17, 
WRITE!NOUT, 530 XMCM!6), YMQM!6) 
F0RMAT!'+* ,83>,*KURT0SIS* ,1X1PE14. 



CALL 

CALL 



REPEAT 
OUTPUT!CELL( 1,18) 



6 ,3XE14.6) 
YMAX -17. 



* DELTA) 



scale, 18, 

WR ITE!N0UT,54C) X!N), YMAX 
540 FORMAT! * + * , 83 X , * MAX IMUM * , 2X1PE14. 6 , 2 XE 14. 6 ) 

CALL REPEAT 

WRITE!N0UT,55C ) X!l), YMIN 
550 FORMAT! * + * ,83X, * MIM IMUM *,2X1PE14.6,2XE14.6) 

Z = YMAX - 17. * DELTA 
DC 560 J = 19,32 
2 = Z - DELTA 

CALL 0UTPUT!CELL! 1, J) , SCALE, J, 2) 

CALL REPEAT 
560 CONTINUE 

WPITE!N0UT,31O 
XLAEEL(l) = X!l) + OELTAX 
FOUPO = * DELTAX 

CO 570 1=3,8 

XLAeELU) = XLABEL(I-l) + FCURD 
570 CONTINUE 

WRITE!N0UT,58O!XLAeEL( I ) , I =1 , 7, 2 ) , ! XLA8FL ! J ) , J = 2 , 8 , 2 ) 
580 FORMAT! 16X,8! * I * ,7X)/ 12X ,4! IPEIO .3 , * I *)/ 20X,4!5 
C10.3,6X) ) 

KEY!1) = 0 
CO 590 1=2,16 
KEY !I ) = SCALE * ^ - 1 
590 CONTINUE 

WRITE !NOUT, 600 ! K EY ! I ) ,I = 1 , 1 6 ) 

600 FORMAT!//50X, 'KEY* , //* SYMBOL PRINTED*, 12X 
1- = = = = T > 

2 F H F F*/ * +* ,38X, * . 

3+X<XVE E 

4/ *+*,56X,*. . - - ■ 

5/ T X*/ * +* ,nox, *T*// 



I — 
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6 • NO. OBSERVATIONS’ , 15, 1516) 

RETURN 

END 

FUNCTION WEIGHTd, K) 

INTEGER WEIGHT 

THIS FUNCTION SUBPROGRAM COMPUTES Tf-E THE WEIGHT FUNCT 
ION FOR THE I-TH RANK IN A SIEGEL-TUKEY TEST WITH A SA 

-mple of N. 

K = N00(I,2) 

WEIGHT = I + I - K 

IF(I .GT. N) WEIGHT = 4 * N - WEIGHT + 2 * (1 - K) 

RETURN 

END 

SUBROUTINE OUTPUT (COUNTS, SCALE, L, Y) 

THIS SUBROUTINE OUTPUTS A SINGLE LINE OF THE BIVARIATE 
HISTOGRAM ON THE LINE PRINTER 
ON THE LINE PRINTER 

INTEGER*2 COUNTS, SCALE, ZERO, ONE 
LOGICAL*! LINE, CHAR 

DIMENSION LINE(64,4), CHAR(4, 16), COUNTS(32), N(16), 
CCHA (16) 

EQUIVALENCE (CHAR(1,1), CHA(D) 

DATA ONE/1/, ZSRO/0/, NQUT/6/, N/3* 1 , 3* 2 , 2*3, 2 , 6 

C*3, 4/ 

DATA CHA /• ’,’=1 

1 *,'=*. ’, •=+. ','TX ',•><- ’,’TX- 

2’,’AVT ',’HE= ','HE/ ', 'HET ’,’HEXT'/ 

LINES = 1 
K = 1 

CO 20 J = 1 ,32 

INDEX = CGUNTS(J) / SCALE + ONE 

IFdNDEX .EQ. 1 .AND. CCUNTS(J) .GT. ZERO) INDEX = 2 
IF(NdNDEX) .GT. LINES) LINES = NdNCEX) 

DO 1C 1=1,4 

LINE(K,I) = CHAR( I, INDEX) 

LIN£(K+1,I) = CHAR ( I ,IN0EX) 

10 CONTINUE 
K = K + 2 
20 CONTINUE 

IF (V0D(L-2,4) .EQ. 0) GO TO 40 
ENTRY REPEAT 

WRITE(N0UT,30) ((LINE(J,I), J=l,64), 1=1, LINES) 

30 FCRNAT(13X, • 1 ',64A1,’| •/(• + ', 13X , 64A 1 ) ) 

GO TO 60 

40 WRITE(NOUT,50) Y, ((LINE(J,I), J=l,64), 1=1, LINES) 

50 FOPMAT( 1X1PE10.3, ' *',64A1,’*'/ ( • 4 • , 13X , 64A 1 ) ) 

60 RETURN 
END 



SUBROUTINE MOMENT (WORK, XMOM) 

PURPOSE J 

FINO'thE moments of the data ARRAY X WHEN NOT ALL 
OF the ARRAY IS AVAILABLE AT ANY GIVEN TIME. 

USAGE: 

DATA POINTS ARE PASSED TO THE SUEPOUTIME "M" AT A T 
-IME WITH THE RUNNING TQTALS KEPT IN the DOUBLE PRECIS 
-ION WORK AREA VECTOR ’WORK*. PRIOR TO THE FIRST CALL 
,W0RK(1) SHOULD BE SET TO 0.0. DATA l$ THEN PASSED BY 
THE CALL 

CALL ACCUM(X, M, WORK) 

AN INITIAL mean ESTIMATE IS KEPT IN W0RK(4) AFTER THE 
FIRST 7 DATA POINTS HAVE BEEN PASSED. HIGHER CENTRAL 
MOMENTS ARE FOUND USING THIS ESTIMATE. ONCE AT LEAST 7 
X VALUES HAVE BEEN ACCUMULATED, RESULTS MAY RE OBTAINS 
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-0 FROM 

call MCMENT( work , XMQM) 

PESLLTS ARE RETURNED AS FOLLOWS: 

>N0M<1) the sample mean 

>M0M(2) the sample variance (UNEIASED ESTIMATOR) 
XM0M(3) THE SAMPLE THIRD CENTRAL MOMENT (LNEIASED 
ESTIMATOR) 

>MOM(4) THE sample COEFFICIENT CF SKEWNESS 
XM0M(5) THE SAMPLE FOURTH CENTRAL MOMENT (UNBIASED 
ESTIMATOR) 

>MCM(6) THE SAMPLE COEFFICIENT CF KURTQSI S 



IMPLICIT REAL*8 (D) 
pcAi.3^8 WORK 

DIMENSION W0RK(8)» XMCM(6), X(M) 

C 

IF(^ORKd) .LT. 7.000) RETURN 

C FINC CORRECTION FACTORS 

CMEfiN = (WQRK(2) + WORK (3 ) ) / WORK(l) 

DC = OMEAN - W0RK(4) 

CC2 = DC =* DC * WORK(l) 

C DETERMINE CORRECT CENTRAL MOMENTS 

DM2 = W0RK(5) - 0C2 

CM3 = W0RK(6) + W0RK(7) - DC * (3. CO W0RK(5) - 0C2 - 
CCC2 ) 

0M4 = W0RK(8) -DC * (4. DO * (W0RK(6) + W0RK(7)) - DC 
1 *(6.C0 ♦ WCRK(5) - 3. DO * CC2)) 

C CORRECT ESTIMATORS IN WORK ARRAY 

W0RK(4) = CMEAN 
W0RK(5) = DM2 
W0RK(6) = 0. 

WCRK(7) = 0. 

IF(CM3 .GT. 0.) W0RK(6) = DM3 
IF(CM3 .LT. 0.) W0RK(7) = DM3 
W0RK(8) = 0M4 

C RETLRN STATISTICS 

XMCM(l) = OMEAN 
CFAC = WORK( I ) - 1 .DO 
XMCM(2) = DM2 / OFAC 
CFAC = OFAC * (WORK(l) - 2. DO) 

XMCM(3) = WORK(l) ♦ DM3 / OFAC 
XM0M(4) = XM0M(3) / (XM0M(2) ♦ 

CFAC = DFAC * (WORK(l) - 3.00) 

XM0M(5) = (((WORK(l) - 2.00) * 

1 - ( 6,00 * WORKd ) 

2WQRK(l)) /DFAC 

XM0M(6) = XM0M(5) / (XMCM(2) * 

RETURN 
C 
C 

ENTRY ACCUM(Xt M, WCRK) 

C 

IF(UORKd) .GT. 6.00) GO TO 210 

C ACCUMULATE FIRST 7 DATA POINTS 

N = WORK( 1) + 2. DO 
WORK(l) = WORKd) + M 
NP = N + M - I 
IF(NP ,GT. 8) NP = 8 
J = 1 

CO 201 I=N» NF 
WORK( I) = X( J ) 

J = J + 1 

201 CONTINUE 
IF(NP .LT. 8) RETURN 

C CBTAIN INITIAL MEAN ESTIMATE 

CSUMP = 0. 

CSUMM = 0. 

CO 202 I = 2t8 

I=(WORK(!) .LT. 0.) DSUMM = DSUMM + WORK(I) 

IF(^^aRK(I) .GT. 0.) DSUMP = DSUMP + UORK(I) 

202 CONTINUE 



SQRT(XMOM( 2) ) ) 

WORKd) + 3.^0) * 0M4 
- S.CC) * DM2 ^ CM2/ 

XM0M(2)) - 3.0 
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IF(J .GT. M) GO TO 20A 
CO 203 I=J,N 

IF(X(I) ,LT. C.) OSUMM = DSUMM + X(I) 
IF(X(I) ,GT. C.) OSUMP = OSUMP + X(I) 
2C3 CONTINUE 

204 CME4N = (CSUMN + DSUMP) / WORK(l) 

C-^ SET yp WORK AREA 

D2 = C. 

C3M = 0. 

D3P = 0. 

C4 = 0. 

CO 205 I=2t8 

CIFF = WORKd ) - OMEAN 

CF2 = DIFF ♦ CIFF 

C2 = 02 + DF2 

IP(CIFF .LT. C.) D3N = C3M + DIFF ♦ CF2 

IF(CIFF ,GT. C.) D3P = C3P + DIFF ♦ DF2 

C4 = 04 + D«=2 * OP2 

205 CONTINUE 

IF(J .GT. H) GO TO 207 
CO 206 I=J,M 
CIFF = X( I) - CMEAN 
CP2 = DIFF ♦ CIFF 
C2 * D2 + DF2 

IF(CIFF .LT. 0.) 03N = C3M + DIFF * CF2 

IF(CIFF .GT. C.) D3P = C3P + DIFF * CF2 

C4 = 04 + DF2 ♦ 0*=2 



206 CONTINUE 



207 V»0RK(2) 


S 


DSUMP 


WORK (3) 


= 


CSUMM 


WOP K(4) 


S 


DME AN 


W0RK(5) 


= 


02 


WCRK(6) 


= 


03P 


W0RK(7) 


= 


D3M 


W0RK(8) 


s 


D4 


RETURN 







C ACCIMLLATE DATA 

210 WCRK(l) = WORKd) + M 
CO 211 I=1»M 

IF(X(I) .GT. C.) W0RK(2) = W0RK(2) + X(I) 

IF(X(I) .LT. 0.) WORKO) = WORK(3) + X(!) 

CIFF = X( I) - W0RK(4) 

CF2 = DIFF ♦ DIFF 
W0RK(5) = W0RK(5) + OF? 

IF(CIFF .GT. C.) W0PK(6) = W0RK(6) 4 DIFF ♦ OF2 

1F(CIFF .LT. C.) WGRK(7) = WORK(7) + DIFF « DF2 

W0RK(8) = W0RK(8) + CF2 * 0F2 

211 CONTINUE 
RETURN 
END 

C 

C 

SUBPCUTINE SCRTON (ON, II, JJ, WITF) 

THIS SUBROUTINE SORTS THE VECTOR ON INTO 
ASCENDING ORDER FROP CN(II) TO ON(JJ), CARRYING 
ALONG THE CORRESPONDING ELEMENTS OF THE VECTOR 
WITH. 

INTEGER HI, HISTK, FIKEEP, STKPT, WITH 
DIMENSION QN(JJ), WITH(JJ), HISTK(16), L0STK(16) 

STKPT = I 
LOKEEP = II 
FIKEEP = JJ 

10 IF (LOKEEP .GE. HIKEEP) GO TO 100 
20 LO = LOKEEP - 1 

M = HIKEEP + 1 

MIDDLE = (LOKEEP + HIKEEP) / 2 
IL = LOKEEP 

IH = HIKEEP 

IF ( ON(MIDDLE) .GT. ON(HIKEEP) ) IF = MIDDLE 
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IF 


( 


ON(MIOOLE) 


.LT. ON(LOKEEP) 


) 


U = MIDDLE 


IF 


( 


0N( HIKEEP) 


.LT. ON(IL) 


) 


IL = HIKEEP 


IF 


( 


QN(LOKEEP) 


.GT. ON(IH) 


) 


IH = LOKEEP 


ISUE 


= LOKEEP + 


MIDCLE + HIKEEP 


— 


IH - IL 



TEST = ON(ISUe) 
ca fQ 50 



30 IF (LO .EQ. HI) GO TO 50 
T£MF = CN(HI) 

QN(M) = CN(LC) 

CN(LO) = TEMP 
ITENP = WI TH( HI ) 

V^ITMHI) = WITH(LO) 

WITMLO) = ITEMP 
50 HI = HI - 1 

IF ( ON(HI) ,GT. TEST) GO TO 50 
70 LC = LO + 1 

IF ( CNILO) .LT. TEST) GO TO 70 

IF ( LO .LE. HI ) GO TO 30 

!F(HI - LCKEEP .LE. HIKEEP - LO) GO TO 90 

LQSIK(STKPT) = LOKEEP 

HISIK(STKPT) = HI 

STKPT = STKPT + 1 

LOKEEP = LO 

GO TO 110 

90 LOSTK(STKPT) = LO 

HISTK(STKPT) = HIKEEP 
STKPT = STKPT + 1 
HIKEEP = HI 
GO TO 110 

100 STKPT = STKPT - 1 

IF( STKPT .EQ. 0) RETURN 
LCKEEP = LOSTK(STKPT) 

HIKEEP = HISTK(STKPT) 

110 IF ( HIKEEP - LOKEEP .GE. 11) GO TO 20 
IF ( LOKEEP .EQ. II) GO TO 10 
LOKEEP = LOKEEP - 1 
120 LCKEEP = LOKEEP + I 

IF ( LOKEEP .EQ. HIKEEP ) GO TO ICO 

IF I ON(LCKEEP) .LE. CN(L0KE£P+1) ) GO TO 120 

TEMP = CN(LOKEEP + 1) 

LQ = LCKEEP 
130 CN(L0+1) = CN(LO) 

LC = LO - 1 

IF ( TEMP .LT. CN(LO) GO TO 130 
CN(L0+1) = TEMP 
IW = LCKEEP + 1 
N = LCKEEP - LO 
ITEMP = WITH( IW) 

CO 140 IMOVE = !♦ N 
mTH(lW) = WITH(IW-l) 

IW = IW - 1 
140 CONTINUE 

V»ITH(IW) = ITEMP 

60 TO 120 

END 
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FCP GFN5R4T1NG EIVARUTE RANDOM V ECTORS < Y, Z ) , CA LL 
FCLLCWING SUBROLTINE. 

1. LNIFORM(0,U MARGINAL DISTRIBUTION. 

CALL LEWUFN(CR,NV.Y,ZtIS); POSITIVE CORRELATION BY 
LAWERENCE ftND LEV.IS METHOD. 

CALL GAVUFMCR ,NV,Y, Z, IS ) ; NEGATIVE CORRELATION BV 
TRANSFORMATION OF GAVER’S EXPONEMTIAL. 

2. EXPCNENTIAL(MEAN=1.0 ) MARGINAL D I S TR I BUT 1 0 M . 

CALL MAREXF(CR,NV,Y,Z.IS) ; POSITIVE CORRELATION BY 
MARSHAL ANC OLKINS METHOD. 

CALL GAVEXFICR ,NV,Y,Z, IS ) ; NEGATIVE CORRELATION BV 
GAVER’S method. 

WHERE 

CR; GIVEN CCRRELATICN. 

NV; NLMBEP OF RANDOM VECTORS BE GENERATING. 

Y; THE first ELEMENT OF VECTOR WITH CIMENSICN NV. 

2; THE SECOND ELEMENT OF VECTOR WITH DIMENSION NV. 

IS; INITIAL RANDOM SEED. 



SUBROUTINE LEWUFM ( CR ,NV »Y , Z , lU ) 

DIMENSION Y (NV) ,Z(NV) ,U( 3) 

IF (CR.LT .0.0) GO TC 20 

A = SCRT(36.0+( 12.0=i'CR ) + ( 16.0*(CR=i'*2 .C ) ))-6 .0 
EETA=A/(2.0*CP) 

ALPf = 2.0-a.0/8ETA) 

F=( l.C-BETA) /(1.0-BETA+(ALPA»BETA ) ) 

R = 1 .0-ALPA 
DC ICO 1=1, NV 
CALL RANDOM! IU,U,2) 

V( I ) = U(1 ) 

IF( L(2) .LE.R) GO TC 50 
W= (U(2)-R )/ ( 1.0-R ) 

IF(W.LE.P) go to iO 

Z(I) = (Y(I)**BETA)*((W-P)/(1.0-P))^’*(P*EETA) 
GC TC 100 

10 2(1 ) = (Y(I)**6ETA)*(W/P) 

GO TO 100 
5C V=U(2)/R 

IF(V.LE.P) GO TO 60 

Z( !) = (( W-P)/( 1.0-P) )^*(R*BETA) 

GC TC 100 
6C Z(I)=W/P 
100 CONTINUE 
PETLRN 

20 WRn£(6,30)CR 

30 FCRMAT( 5X, 'THIS METHCO NOT ALLOW CR=',F7.3) 
fftLRN 
ENC 



SUBROUTINE GAVUFM(CP,NV ,Y,Z,I S) 

C 

DIMENSION Y(NV),Z(NV),WK(1),G(100) 

DOUELE PRECISION DXS 

CXS=IS*0.8 

IU=IS*0.5+123A5 

IF (CP.GT .0.0 ) GO TO 20 

P=-2.0*CR 

C=1 .C-P 

DC 100 1=1, NV 

CALL GGEOTICXS ,1 ,0 ,WK1 , IG) 

CALL RANDOM! IL,U, 1 ) 

PN=1.0/IG 

V(I ) = (0=»(U=^'*RN) )/(1.0-(P*(U*«RN) ) ) 
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CALL RANDCMdL'tG, IG ) 

V' = 1.0 

CG 1C K=1,IG 
K = W=*G( K) 

1C CGMINLE 

2( I )=W»«( 1.0-P) 

ICO CONTINUE 
RETURN 

20 WRITE(6,30)CR 

30 FCRNAT(5Xt*THIS METI-CD NOT ALLOW CR=*,FT.3) 
RETURN 
END 
C 
C 

c 

c 

SUBROUTINE M AR EXP ( CR , NV , Y, Z , I X ) 

C 

DIMENSION Y(NV),Z(NV),E(3) 

IP (CR.LT.0.0) GO TO 20 
R12=(2.C*CR)/(CR+1.C) 

Rl=1.0-R12 
R2 = P1 

DO 5C I=1»NV 
CALL EXP0N(IX,E,3) 

Ed ) = E( 1) /R1 
E(2 ) = E( 2) /R2 
£(3 )=E(3)/R12 
V( I)=AMIN1(E(1) ,E(3)) 

Z( I )=AMIN1{ E( 2) »E( 3 n 
50 CONTINUE 
RETURN 

20 WRITE(6,3C)CR 

30 FCRNAT(5X,'ThIS METHCD NOT ALLOW CR=’,F7.3) 
RETURN 
END 



C 

C 

SUBROUTINE GAV EXP ( CR , NV » Y, Z , I S ) 

C 

DIMENSION Y(NV) ,Z( NV) ,WK1( 1) , WK2( ICC) 

CCUBLE PRECISION DS E EDI ,0SEED2 » OS E E C 3 
CSEED1=IS«0.2+123A56 
DSEEC2=IS*0.5d234 
CSEED3 = IS#0 .9 + 123 

IF<CR.LT.-0. 5 .OR. CR.GT.0.0) GO TC 20 
C = -2 .0*CR 
P=1.C-Q 
CO ICC 1=1, NV 

CALL GGEOT( CSEEDl ,1,P,WK1, IG) 

A=IG 
E = P 

CALL GGAMS(CSEE02,A,e, 1 ,WK2,GA) 

Z( I )=GA 

CALL GGUBS( CSEED3, 1,U) 

RN = 1.0/ IG 

XI )=ALOG (( 1. C/(D*U**RN) )-( C/F) ) 

100 CONTINUE 
RETURN 

20 VRITE(6,3C)CR 

30 FGRWAT(5X,»THIS METFCD NOT ALLOW CR=',F7.3) 
RETURN 
END 
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